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 ii 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. iii 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 iv 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 v 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 33 LU EI ∞ =Π μ , cf=Π5 , D χ=Π 6 as Variables... 95 vi 4.3.2 Subset of Space State with 33 LU EI ∞ =Π μ , cf=Π 5 , θ=Π 7 as Variables ... 96 4.3.3 Subset of Space State with μ ρ DU∞=Π1 , 33 LU EI ∞ =Π μ , and D χ=Π 6 as Variables ................................................................................................................... 98 4.3.4 Subset of Space State with μ ρ DU∞=Π1 , 33 LU EI ∞ =Π μ , and 7 θΠ = as Variables ................................................................................................................... 99 4.3.4 Subset of Space State with, D L=Π 2 , 33 LU EI ∞ =Π μ , and D χ=Π 6 as Variables ................................................................................................................................ 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 vii 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 viii 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. smcfm 100508.0 =× . ..................................................... 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 10=DL , o40=α at 40Re = ...................................................................................... 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 ix Table 5.3 Maximum and RMS difference between simulations and curve fit. .............. 126 Table 5.4 Comparison of LC and DC between simulations and curve fits for a specific simulation case that had not been used to derive the curve fit ( o45=α , 7=DL , 7Re = ) ........................................................................................................................... 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 x 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 xi 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: Re−PC 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 mmr 2.0= for Integra T fabric. ........................................................................................................................................... 49 xii Figure 2.20: Spatial variation of ZD velocity when averaged over an ellipse with semi- major axes mm)6.0,5.0,4.0,2.0(= and semi-minor axis mm02.0= . 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 semi- major axes mm)6.0,5.0,4.0,2.0(= and semi-minor axis mm02.0= . 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 xiii Figure 3.10: Re−PC 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: Re−PC 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 xiv 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. 11 =Π , 402 =Π , 25.13 =Π , 84 =Π , 5 0.2Π = , o57 =Π ................. 94 Figure 4.16: Fiber hang-up on a cylinder for different values of flexibility with 2.0=cf , and 1Re = . The fiber stiffness in figures (a), (b), and (c) is respectively [ ] 212 .105,25,80 mNEI −×= . .............................................................................................. 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 2.0=cf , and 1Re = . The fiber stiffness from left to right is [ ] 212 .102,40,140 mNEI −×= , and ]40,16,13[ ooo=θ .............................................................................................................. 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 [ ] 212 .102,80 mNEI −×= from left to right. .................. 98 Figure 4.21: The hung-up region as a function of 1Π , 3Π , and 6Π . .............................. 99 xv 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 40Re = , showing the negligible influence of the wall........................................................................................................ 109 Figure 5.3: Velocity profile along the dashed line at 1Re = , 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 2=DL , 40Re = , o90=α . Note the highly expanded ordinate axis........................................................... 112 Figure 5.7: Drag coefficient versus number of mesh volumes for 5=DL , 1Re = , o10=α . Note the highly expanded ordinate axis........................................................... 112 Figure 5.8: Drag coefficient versus Reynolds number ( 40Re0 ≤p ) for a cylinder normal to the flow with DL as a parameter. ................................................................ 114 Figure 5.9: Comparison of drag coefficient between simulations and Cox, 1Re = . ..... 115 Figure 5.10: Streamlines around a cylinder with 20=DL at different inclination angles. The Reynolds number is 40. ........................................................................................... 117 Figure 5.11: Streamlines around a cylinder with 2=DL 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 xvi Figure 5.13: Nondimensional pressure distribution along a cylinder with 20=DL for different inclination angles at 40Re = . (a) Front pressure (b) Rear pressure. ............... 120 Figure 5.14: Forces on a finite cylinder with ( =DL 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) 1Re = , (b) 40Re = . ................................................... 122 Figure 5.16: Nondimensionalized torque as a function of the inclination angle for different aspect ratios. (a) 1Re = , (b) 40Re = . ............................................................. 123 Figure 5.17: Non-dimensional drag coefficient ( =DL 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 ( =DL 2, 5, 10, 20) versus the angle of inclination with Re as a parameter................................................................... 125 xvii 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. xviii 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 xix 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. xx 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. xxi 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. xxii 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. 1 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 mμ150 in diameter, with a pitch of mμ400300 − . Typical pulp fibers have a diameter of mμ4020 − and are mμ3000100 − long, and pulp fines are fiber fragments less than mμ100 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 twin- wire/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 2 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 paper- side 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. 3 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 single- layer 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.” 4 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] Investigated the effect of long knuckle orientation on drainage and sheet formation. (single-layer fabric) Faster drainage and better release when fibers oriented across the long knuckles. Beran (1979) [15] Introduced a fiber support index (FSI) based on mean fiber length and the average number of supports per fiber. (single-layer fabric) FSI can determine optimum fabric designs for a particular paper machine. Higher FSI corresponds with papers with finer grades. Johnson (1984) [8] Introduced a drainage index (DI) based on fabric air permeability, CD mesh counts and the FSI. DI is used to assess drainage performance. Higher values of DI reflect better drainage. Johnson (1986) [9] Using DI predicted retention and drainage behaviour of multilayer fabrics. Inverse relationship of retention to MD frame length; proportional relationship of drainage to DI. Adanur (1994) [5] Investigated effects of fabric characteristics on sheet formation and properties. An increase in DI and FSI resulted in an increase in sheet strength; An increase in mesh count reduced sheet thickness. Helmer et al. (2006) [16] Introduced formation index (FI) as a function of “crowding number”, “Froude number” and “aspect ratio of flow depth to fiber length”. Characterized dewatering as “early” and “late” drainage velocity; Higher late drainage rate refers to higher FI and higher topside anisotropy. Xu et al. (2009) [17] Studied the effect of internal structure of triple fabrics on drainage behaviour. An increase in the center plane resistance increases the total resistance only at lower basis weight. 5 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 Re−K . This coefficient, K , is usually defined as the ratio of static pressure drop across the screen against the dynamic pressure of the upstream ( 2 2 1 ∞ Δ= U pK ρ ). It should be noted that researchers have used different definitions of Reynolds number for screens: (i) Based on the filament diameter as μρ dUd ∞=Re . (ii) Based on the ratio of projected open area of the screen to the total cross- sectional area, “porosity” β , as βdp ReRe = . (iii) Based on the surface area exposed to the flow per unit volume, S , ( )Sdds ReRe = . (iv) Based on the wetted perimeter of the orifice where the flow is most constricted, W , and the filaments’ center-to-center distance, l , ( ) dw Wdl Re4Re 2= . 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. 6 Table 1.2 Screen/grids resistance-flow rate relationships from experiments Author Description Robertson (1950) [23] By making an analogy to the flow through orifices and nozzles, presented the resistance, discharge coefficient, over a range of 4000Re20 << d . In contrast to a single orifice, at high flow rates the discharge coefficients could be greater than unity. Wieghardt (1953) [24] By means of experimental measurements over the range 6000Re60 << s , he correlated the pressure drop coefficient, K , to sRe and β . Ingmanson et al. (1961) [25] Modeled the flow through screens with the flow through porous media. Using experimental values, they proposed a generalized equation in the form Css BAf ReRe 1+= − for 1500Re3 << s . Screen structural properties are considered in the constants of the correlation. Armour and Cannon (1968) [26] Modeled screens as a thin packed bed and proposed the correlation BAf e += −1Re for 1000Re1.0 << e , applicable to diverse groups of screens. Pedersen (1969) [27] Compared the flow with orifice flow and presented a two-variable discharge coefficient analytically, ( )Wd fC Re= . He neglected the flow interaction between neighboring orifices and considered cylindrical filaments in both the Machine Direction (MD) and Cross Machine Direction (CMD). Wakeland and Keolian (2003) [28] Measured the low-frequency oscillating flow across square-mesh screens in the range 400Re002.0 << d and correlated K , to dRe and β . 7 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 BWd AC Re= . Wang et al. (2007) [35] Following the work of Pedersen derived a discharge coefficient, dC , for filaments of elliptical cross-section ( )baar = . They provided ( ) DrWBrd CaLnAaC += Re . For the same Re, increasing ra decreased dC . Huang (2003) [36] He simplified fabrics as rows of dissimilar cylinders studied the 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. 8 Table 1.4 Slender particle in a creeping flow of a Newtonian fluid Author Description Jefferey (1922) [38] Found the motion of a neutrally buoyant rigid elliptical particle in a simple shear flow ( )0,0,yU γ&v =∞ and showed the dimensionless period of rotation, γ&T , depends merely on the aspect ratio of ellipse ( )PP rrT 12 += πγ& . Forgacs and Mason (1958) [39] Using the Euler’s equation they showed that the minimum value at which rod-like particles of axis ratio pr and bending modulus bE , may be expected to buckle under shear-induced compress is given by ( ) ( ) 4275.12 PPbcrit rrLnE −≈μγ& . Bretheron (1962) [40] Showed Jeffrey’s analysis to be valid for any axisymmetric body with fore-aft symmetry, provided that an equivalent aspect ratio er is used in place of the actual spheroid aspect ratio, pr . Tillet (1970) [41] 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. Batchelor (1970) [42] Modified the Tillett’s work for non-circular shape of the cross-section of the body. Cox (1970) [43] Cox’s theory allowed also a curvature of the body comparing to the work of Tillet and Batchelor. Khayat and Cox(1989) [44] 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 (2005, 2006) [46, 47] Studied the inertia effect in more detail and found a gentle ( )ReO drift in closed Jeffrey orbits. 9 Table 1.5 Main methods to simulate flexible fibers in a flow field Author Description Hinch (1976) [48] Utilized the slender-body theory of Cox to formulate a flexible inextensible isolated thread in a shearing flow. Gradon et al. (1988-1992) [49,50,51] 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. Yamamoto and Matsuoka (1993-1997) [52-55] 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 Kuhn (1997) [56] Presented a method to describe the dynamics of a flexible body based on the equations of motion for a beam in vibration. They examined a flexible fiber motion in a channel flow with a slot. Klingenberg et al. (1997-2002) [57-60] To eliminate the need for iterative constraint to maintain fiber connectivity in Yamamoto’s method, they proposed a particle-level 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 Green (1998) [62] They applied the Immersed Boundary Method (IBM), based on a mixture of Eulerian and Lagrangian variables [61], to fiber motion in a shear flow. Due to the two-way coupling between phases, they could observe disturbance in the flow caused by the fiber motion. 10 Table 1.5 Main methods to simulate flexible particle in a flow field (cont.) Zhu and Peskin (2002,2003) [63,64] Zhu used the IBM to simulate flapping filaments in a flowing soap film. His results showed the amplitude of flapping increased as the mass of the filament increased, needed for a sustained flapping. In the case of two flapping fibers he found two distinct modes of sustained 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 Shelley (2004) [66] Based on a non-local slender body theory, they introduced a method to simulate interactions of flexible fibers in Stokes flow. They used the forces described by Euler-Bernoulli elasticity theory in lieu of the slender body integral equations. Qi (2006) [68] 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 (2008) [69] In this recent work, Lindstorm extended Schmid’s approach to include the two-way coupling between the fiber and the fluid phase and particle 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. 11 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 sm . The flow around them is approximately that around a circular cylinder of aspect ratio 1002 − with arbitrary yaw angle, and with a Reynolds number, μρ DU∞=Re ( =ρ 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 40Re ≈ , and thereafter becomes unsteady and laminar ( 200Re40 << ) and then for 300Re > 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 4104× [82]. He found that for aspect ratios in the range 40
, unsteady simulations were also done. These simulations produced results within 2% of the steady computations. 40 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 41 coarse grid ( 2h ) to that of the fine grid ( 1h ). As integral quantities are used for error estimations, these representative cells are chosen as: ( ) 3 1 1 1 ⎥⎦ ⎤⎢⎣ ⎡ Δ= ∑ = N i iVN h where iVΔ is the volume of the thi 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 ( 57.2 ′′ ) 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. 42 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. smcfm 100508.0 =× . 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. 43 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 44 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 45 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 (= mm15.0 for fabric 1 and mm125.0 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, 2R , is very high, indicating the quality of the function fit. Furthermore, for both fabrics PC is approximately proportional to 4.0Re− . Interestingly, the exponent on the Reynolds number is almost exactly half way between the creeping flow limit ( 1Re−∝PC ) and the fully developed turbulence limit ( 0Re∝PC ). Since 2VCp P∝Δ and 4.0Re−∝PC and V∝Re , therefore 58Vp ∝Δ . 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 drive 616 125 250440 8 5 =⎟⎠ ⎞⎜⎝ ⎛× cfm through the same fabric. 46 Figure 2.17: Re−PC 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 mmr 4.0> . In broad terms, this means that a disk of radius greater than mm4.0 would be pulled onto the fabric at about an equal rate, irrespective of the lateral location of the disk. 47 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 mmr 2.0< , 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 semi- minor axis of 0.02mm (i.e., roughly corresponding with (0.4, 0.8, 1, 1.6) mm long fibers mμ40 in diameter). The results of the averaging are shown in Figure 2.20 for o0=α (aligned with MD) and in Figure 2.21 for o90=α (aligned with CMD). The contours of o0=α , 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 o90=α 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. 48 (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. 49 Figure 2.19: Spatial variation of circular averaging over mmr 2.0= for Integra T fabric. 50 (a) semi-major axis mm2.0= (b) semi-major axis mm4.0= (c) semi-major axis mm5.0= (d) semi-major axis mm6.0= Figure 2.20: Spatial variation of ZD velocity when averaged over an ellipse with semi- major axes mm)6.0,5.0,4.0,2.0(= and semi-minor axis mm02.0= . The ellipse is aligned with MD. Velocities in meters/second. 51 (a) semi-major axis mm2.0= (b) semi-major axis mm4.0= (c) semi-major axis mm5.0= (d) semi-major axis mm6.0= Figure 2.21: Spatial variation of ZD velocity when averaged over an ellipse with semi- major axes mm)6.0,5.0,4.0,2.0(= and semi-minor axis mm02.0= . The ellipse is aligned with CMD. 52 (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. 53 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. 54 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 Three- Dimensional 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. 55 [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);364- 368, 1998. 56 Chapter 3 The Influence of Machine-Side Filaments and Jet-to-Wire Speed Ratio on the Flow Through a Forming Fabric2 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. 57 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. 58 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. jet velocity fabric velocity =FV =jV α impingement angle Oblique dewatering Perpendicular dewateringFree-jet zone free surface Figure 3.1: Schematic of fabric and impinging jet. 59 In papermaking, the incoming jet impinges on the fabric with angles in the range oo 103 << α , as observed in the frame of reference of the papermachine. In terms of the fabric reference frame, this impingement angle, θ , is: 1cos sintan −× ×= α αθ JTW JTW in which W j V VJTW = is the jet-to-wire speed ratio, which is normally between 0.990 and 1.035 [16]. Therefore, for an impingement angle of o5=α with a 035.1=JTW , the corresponding angle in the fabric reference frame would be o70≈θ . 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. 60 Table 3.1 General specification of Fabric A used in simulations. Mesh/knock (1/inch) Yarn Size (mm) Paper side Machine side Total Paper side Machine side Air Permeability (cfm) MD CD MD CD 146×186 73×93 73×62 440 0.13 0.13 0.21 0.20 Figure 3.2: Hybrid computational mesh scheme. 61 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 62 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 ( 2h ) to that of the fine grid ( 1h ). As integral quantities are used for error estimations, these representative cells are chosen as: ( ) 3 1 1 1 ⎥⎦ ⎤⎢⎣ ⎡ Δ= ∑ = N i iVN h where iVΔ is the volume of the thi 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 63 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 64 (a) paper-side layer (b) machine-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. 65 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. 66 3.3.3 Averaged Velocity Contours Typical pulp fibers are mμ3000300 − 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 paper- side 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. 67 % % 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 (= mm125.0 ), and ν is the kinematic viscosity of the fluid. The Re−PC 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. 68 Figure 3.10: Re−PC 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 paperCP , 27 21 18 17 15 14 machineCP , 42 33 29 27 24 22 PC∑ 69 54 47 44 39 36 bothCP , 64 52 44 42 37 34 % difference 7 2 2 2 5 3 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 oo 9058 ≤≤ θ , corresponding to impingement angles oo 153 <≤ α and 035.199.0 ≤≤ JTW . 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. 69 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 ( 2 2 1 VPCP ρΔ= , where V is the bulk velocity). The pressure coefficient is plotted 70 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: Re−PC curve for Fabric A with impingement angle as a parameter. 71 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. 72 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. 73 [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. 74 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) fibre mat Machine Direction 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. 75 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 er was used in place of the actual spheroid aspect ratio, pr [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. 76 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 particle- level 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 77 dimensionless groups. By creating simulations that involve different sets of dimensionless groups we are able to draw a demarcation line between regions in the 7- dimensional 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 1Re << [22], but for 1Re ≥ , there is no analytical work that describes the flow in detail. Therefore, in anticipation of future studies where 1Re >> , 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. 78 upL downL 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] 275.089.0 Re1064.31 Re0004.0 Re 96.1 Re 8.618.1 D D DD DC −×+−++≈ We conducted our simulations in a channel to save computational effort. A uniformity of the x-component of velocity was attained, typically for DY 10> , whereas the y- component 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 1.0005.0 << λ . The λ−DC 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 5.0Re = the DC value from numerical data in [23] is 79 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 DLup 10= , DLdown 20= and DW 10= 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]. 80 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 each segment i is located at ir r , measured from a fixed reference frame. Each segment is a prolate spheroid with a major axis a2 and a minor axis b2 , as illustrated in Figure 4.6. The orientation vector of each segment is iP r . 81 ds Y X Z Y X Z ix iy iz 1=i Ni = ir r Figure 4.5: Representation of a fiber as a series of N linked rigid bodies. Y X Z i rr iP rJoint Joint i 1+i a2 b2 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. 82 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 the angular velocity ωr and the translational velocity Ur , the hydrodynamic torque and force are given by [22] ( )jjijhi UUAF −= ∞ (4.1) ( ) ∞∞ +−Ω= jkijkljjijhi EHCT ω (4.2) where ∞∞∞ Ω EU ,, 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 ( )[ ]jiijAjiAij ddYddXaA −+= δπμ6 (4.3) ( )[ ]jiijCjiCij ddYddXaC −+= δπμ 38 (4.4) klijl H ijkl ddYaH επμ 38= (4.5) in which HCCAA YYXYX ,,,, are only functions of eccentricity. The internal resistance torque in the joints, jY is a bending torque given by ( ) PBjjbbj nkY ˆ0θθ −−= (4.6) where aEIkb 2= is the bending constant, 0jθ the equilibrium value and PBn̂ 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 321 MFMMr cjjiji ++= rrr& ω (4.7) 4321 i c i c jjlil NTNFNN ++= rrrω (4.8) The scalar matrices in the above equations are ( )jkkjiij bAAbM ~~ 11 −+−= (4.9) 12 −= AM (4.10) ∞−= jjUAAM 13 (4.11) ( )[ ]lkkjijljjijiijil bAAAdbAdarCN ~~~~3 121 −+−= δ (4.12) [ ]jijj AIdarN −= ~4 1 32 (4.13) 83 3 3 8 1 ar N = (4.14) ( )[ ] 3124 8 1 ~3 ar CEHUAAAUA ar d N jjjjkkijj ij i +Ω++−= ∞∞−∞ (4.15) where ∑ = = N i iAA 1 and baar = . The tilde operator is used to perform a cross product 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 b2 , bkb 3πμ and bkb . 4.2.2.1 Fiber in an Unbounded Field When a fiber is in an unbounded flow, there are no contact forces, i.e. cF r and cT r are zero. So, the governing equations reduce to 31 MMr jiji += ωr r & (4.16) 41 ilil NN =ωr (4.17) The rotational evolution of fiber can be directly calculated from Eq.4.17. By substitution of ωr 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 0),( =yxS , 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 84 S Sni ∇ ∇= r r ˆ , 0ˆ.ˆ =ii nt (4.18) it̂ in̂ Y X segment i 0),( =yxS 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 constraint on the normal velocity of the segment that touches the wall, i.e. 0ˆ. =ii nur . 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: ( ) [ ] [ ] ⎩⎨ ⎧ = ==⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛⇒ ⎪⎭ ⎪⎬ ⎫ =− =++= contactsi Nj B F A NFNN nMFMMnr c i j i c jjjij i c jjijii :1 , 0ˆ.ˆ. 421 321 r r rr rrr& ω ω ω (4.19) The unknowns are now ωr and cFr ; and cFr 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 of friction i.e. nf t FF rr μ= . In contrast, if the fiber staples to the wall the two reaction forces are independent and should only satisfy the inequality ∑∑ < nifti FF rr μ , where 85 the summation operator is performed over all segments in contact with the wall. The procedure solution is as follows: (1) First, assume that the fiber remains stationary. An extra constraint of 0ˆ. =ii tur 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, ∑∑ < nifti FF rr μ . (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. (4-2) The components of contact forcers are related by nf t FF rr μ= (4-3) Find the velocities and forces from the systems of Eq.4.19. (4-4) Check the velocity if the right-sliding condition is satisfied 0ˆ. [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, 200Re60 << , 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:1771- 1781. 106 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 4104× [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 40=DL to 1000=DL , whereas for smaller aspect ratios the onset of laminar vortex shedding was delayed: 100=DL 5.47Re =⇒ c , 30=DL 5.49Re =⇒ c and 20=DL 5.51Re =⇒ c . It is shown that critical Reynolds numbers fit well into the relation 2 3108.14.47Re − ⎟⎠ ⎞⎜⎝ ⎛×+= D L c [11]. Szepessy and Bearman found that aspect ratio also has an effect on the shedding frequency. For example, 19.0=St for 4105.4Re ×< and 7.6=DL reduced to 17.0=St at similar Reynolds numbers but 1=DL [12]. Norberg wrote a review of previous simulations of fluctuating lift and has found that most simulations have covered the range 500Re45 << and 510Re100 << for two- dimensional 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 300Re40 << . The dependence of the critical Reynolds number and the Strouhal number on DL 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 51 ≤≤ DL [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 107 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 310Re160 << . 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 401− . 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 900 ≤≤α for 40Re1 ≤≤ . 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. 108 5.2 Numerical Method 5.2.1 Meshing and Boundary Conditions Cylinders with four different aspect ratios ( =DL 2, 5, 10, and 20), with angles of inclination 900 ≤≤α 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 L10 upstream to L15 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 10=DL , o50=α , 40Re = and in Figure 5.3 for the case 10=DL , o50=α , 1Re = . As one approaches the lateral boundaries of the domain, typically for LY 2> , the velocity falls within 0.1% of the free stream velocity. Note that the velocity profiles illustrate the dominance of the viscous forces for 1Re = , and the inertial forces for 40Re = : near the cylinder the flow is accelerated for 40Re = , whereas it is decelerated for 1Re = . Figure 5.1: Computational domain indicating the inlet, outlet, and slip walls. 109 Figure 5.2: Velocity profile along the dashed line at 40Re = , showing the negligible influence of the wall Velocity extracted along the dashed line. ∞U Figure 5.3: Velocity profile along the dashed line at 1Re = , 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, 110 whereas wedge elements are used far from the cylinder. Figure 5.5 depicts a close-up view of the unstructured mesh around the cylinders. L2 L14 L9 ∞U 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 111 unsteady vortex shedding normally begins when 50Re ≥ for cylinders of aspect ratio 20≤DL , almost all simulations were run under steady flow conditions. We ran some simulations both steady and unsteady for the case 10=DL , o40=α , 40Re = . 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 10=DL , o40=α at 40Re = 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. 112 Figure 5.6: Drag coefficient versus number of mesh volumes for 2=DL , 40Re = , o90=α . Note the highly expanded ordinate axis. Figure 5.7: Drag coefficient versus number of mesh volumes for 5=DL , 1Re = , o10=α . 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 113 for error estimation, it is appropriate to choose average global representative cells 1h , 2h , and 3h as ( ) 3 1 1 1 ⎥⎦ ⎤⎢⎣ ⎡ Δ= ∑ = N i iVN h Where iVΔ is the volume of the thi 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 DL (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. 114 Figure 5.8: Drag coefficient versus Reynolds number ( 40Re0 ≤p ) for a cylinder normal to the flow with DL 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. Re Simulation DC ( ∞→DL ) Ref [25] Difference (%) Ref [24] Difference (%) 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 DC on DL , consistent with the creeping flow theory of Cox [26], in which the total force F r acting upon a circular cylinder of finite length at rest in a uniform flow of velocity U r is given by jiji UA LF ⎟⎠ ⎞⎜⎝ ⎛= 2 μ 115 where ijA is a tensor given by ( 1=i refers to the direction along the cylinder axis) ( ) ( ) ( ) ( ) ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ ≠ == +− == +− = ji ji D L ji D L Aij , 3,2, 2ln 2 12ln 8 1, 2ln 2 32ln 4 0 π π For 1Re = (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 DL to the creeping flow theory of Cox. Figure 5.9: Comparison of drag coefficient between simulations and Cox, 1Re = . 116 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 LDU FC DD 2 2 1 ∞ = ρ , LDU FC LL 2 2 1 ∞ = ρ 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 20=DL at different inclination angles. The flow Reynolds number is 40. ( o90=α ) ( o70=α ) 117 ( o50=α ) ( o20=α ) Figure 5.10: Streamlines around a cylinder with 20=DL at different inclination angles. The Reynolds number is 40. For a cylinder with 20=DL normal to the flow ( o90=α ), 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 D2.2 . Our simulations, at a Reynolds number of 40 and for a finite cylinder with 20=DL , show a separation bubble at the middle of the cylinder with a length of D8.1 . 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 o50≤α , 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). 118 ( o90=α ) ( o40=α ) Figure 5.11: Streamlines around a cylinder with 2=DL at different inclinations angles. The Reynolds number is 40. At lower Reynolds numbers, e.g. 1Re = , 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. ( o90=α , 20=DL ) ( o40=α , 20=DL ) 119 ( o90=α , 2=DL ) ( o40=α , 2=DL ) 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 20=DL , the pressure at the front stagnation line is nearly constant, with 03.1=PC (see Figure 5.13a). By way of comparison, 2D simulations by Sen et al. [28] show a stagnation point pressure coefficient of 16.1=PC at 40Re = . 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 0.48PC = − (see Figure 5.13b), which compares favourably with the values reported by Sen et al. [28] ( 0.48PC = − ) and Baranyi and Lewis [25] ( 0.51PC = − ). As tabulated in Table 2, at a Reynolds number of 40 the asymptotic value of DC , obtained by fitting an exponential curve through a plot of DC versus L D , is 1.55. For this Reynolds number, in 2D simulations Baranyi & Lewis [25] report 1.54DC = and Sen et al. [28] report 1.51DC = . In both cases the agreement may be considered good, in view of the 1% numerical error and the asymptotic function fitting error. 120 (a) (b) Figure 5.13: Nondimensional pressure distribution along a cylinder with 20=DL for different inclination angles at 40Re = . (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 1Re = and 5Re = . Graphs at higher Re are similar, but show a reduced DC . For a given DL and Re , the drag coefficient increases monotonically as the inclination angle increases. Owing to symmetry, the DC versus α graph must have zero slope at both o0=α and o90=α and therefore must have an inflection point somewhere between. In all cases, this inflection point occurs at o45≈α . 121 Figure 5.14: Forces on a finite cylinder with ( =DL 2, 5, 10, 20) versus the angle of inclination with Re as a parameter. a) Drag coefficient b) Lift coefficient 122 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 o45≈α , 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 1Re = . 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) 1Re = , (b) 40Re = . 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). 123 (a) (b) Figure 5.16: Nondimensionalized torque as a function of the inclination angle for different aspect ratios. (a) 1Re = , (b) 40Re = . 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 1Re = and 5Re = . 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 , D D C C ⊥ 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 ⊥,DC on Re and DL is shown in Fig. 8. 124 Figure 5.17: Non-dimensional drag coefficient ( =DL 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 DL CC as a function of α , with DL and Re as parameters. Because the lift coefficient curves show a maximum at angles o45≈α , the ratio of lift to drag coefficient is also a maximum around o45=α . Furthermore, since LC versus α is nearly symmetric with respect to o45=α , whereas DC increases with increasing α , the DL CC curve is asymmetric with respect to o45=α , being higher near o0=α than near o90=α . The asymmetry of the graphs becomes larger as DL rises. 125 Figure 5.18: Lift-to-drag ratio on a finite cylinder with ( =DL 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 40Re = for o10=α or o80=α . At smaller Re and inclination angles oo 8010 << α this discrepancy is limited to less than 10%. The largest discrepancy in drag coefficient occurs for 40Re = and for o30<α . 126 Table 5.3 Maximum and RMS difference between simulations and curve fit. drag coefficient lift coefficient DL 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 o45=α , 7=DL and 7Re = , which had not been used as data for the curve fitting. The values of LC and DC 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 LC and DC between simulations and curve fits for a specific simulation case that had not been used to derive the curve fit ( o45=α , 7=DL , 7Re = ) Simulation Curve fit Difference % LC 0.65 0.63 3.1 DC 2.78 2.93 7.3 According to the Independence Principle, LDUCF NDN 2 2 1 ρ= , where NU is the velocity component normal to the cylinder axis, i.e. sinNU U α∞= . 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 α , DL , and Re . In this table the simulation drag is computed from the best fit curves. The Independence Principle drag is computed with DC based on the Reynolds number of the normal component of velocity ( Re /NU Dρ μ= ). The Independence Principle is accurate to within about 5% for o70≥α , overestimates the drag by 15% or more for o45=α , and can be more than 50% 127 in error for o30≤α . The error of the Independence Principle is greater at larger Reynolds numbers and at larger values of DL . 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. Cylinder Geometry o10=α 30Re = 15=DL 30α = o 30Re = 15=DL 45α = o 7Re = 7=DL 45α = o 30Re = 7=DL 45α = o 30Re = 15=DL 70α = o 2Re = 2=DL 70α = o 30Re = 15=DL Sim IP F F 0.60 0.61 0.87 0.80 0.74 1.00 0.94 Table 5. 6 Lift coefficient curve fit )4sin(Re,)2sin(Re,Re,, 42 ααα ⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛=⎟⎠ ⎞⎜⎝ ⎛ D LA D LA D LCL 9276.11 =λ 9124.12 −=λ 4270.03 =λ 2 1 2 Re Re, δδ +=⎟⎠ ⎞⎜⎝ ⎛ D LA )exp( 3211 D L D L λλλδ −+=⎟⎠ ⎞⎜⎝ ⎛ )exp( 6542 D L D L λλλδ −+=⎟⎠ ⎞⎜⎝ ⎛ 3870.04 =λ 4420.05 −=λ 5050.06 =λ 0330.07 =λ 1366.08 −=λ 3912.09 =λ 43 2 4 ln(Re) δδ += A A )exp( 9873 D L D L λλλδ −+=⎟⎠ ⎞⎜⎝ ⎛ )exp( 1211104 D L D L λλλδ −+=⎟⎠ ⎞⎜⎝ ⎛ 0515.010 =λ 0588.011 −=λ 0382.012 =λ 128 Table 5.7 Drag coefficient curve fit ⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛=⎟⎠ ⎞⎜⎝ ⎛ ⊥ Re,)2cos(Re,Re,, 02 , D LA D LA D L C C D D αα 0414.01 −=γ 07010.02 =γ 9342.03 =γ ( ) 212 RelnRe, ββ +=⎟⎠ ⎞⎜⎝ ⎛ D LA ( )DLDL 3211 exp γγγβ −+=⎟⎠⎞⎜⎝⎛ ( )DLDL 6542 exp γγγβ −+=⎟⎠⎞⎜⎝⎛ 2377.04 −=γ 2685.05 =γ 1697.06 =γ 0470.07 −=γ 0790.08 =γ 8194.09 =γ ( ) 430 RelnRe, ββ +=⎟⎠ ⎞⎜⎝ ⎛ D LA ( )DLDL 9873 exp γγγβ −+=⎟⎠⎞⎜⎝⎛ ( )DLDL 1211104 exp γγγβ −+=⎟⎠⎞⎜⎝⎛ 7609.010 =γ 2577.011 =γ 1583.012 =γ in which 2ReRe, 1, κκ=⎟⎠ ⎞⎜⎝ ⎛ ⊥ D LCD 2 1 1 ηηκ +=⎟⎠ ⎞⎜⎝ ⎛ D LD L 4 3 2 ηηκ +=⎟⎠ ⎞⎜⎝ ⎛ D LD L 7565.181 =η 8586.102 =η 3233.03 −=η 6024.04 −=η 129 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, =DL 2, 5, 10, and 20 oriented to the flow at angles between oo 900 ≤≤α and for 40Re1 ≤≤ 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 DL and Re , the drag coefficient increases monotonically with α , whereas the lift coefficient is approximately parabolic, achieving a maximum near o45=α . The curves of LC and DC against α are similar as one varies DL and Re . Both LC and DC decrease with increasing Re . LC increases rapidly with DL for 10≤DL , and thereafter is nearly independent of DL . DC decreases rapidly with DL for 10≤DL , and thereafter is likewise nearly independent of DL . Curve fits of ( )Re,, DLCL α and ( )Re,, DLCD α 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 o70≥α ), and overestimates the drag force by increasing amounts as α is diminished (error > 15% for o45=α and error > 50% for o30≤α ). 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. 130 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 ( ) ( )100Re1 OO << . 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. 131 [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. 132 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 b p aC Re= , with coefficient of determinations close to 1. 2) At the moderate Reynolds number of the studies, 140Re20 ≤≤ , all calculated powers on the Reynolds number fell in the range 01 ≤≤− b , where 1−=b corresponds to the creeping flow and 0=b to the fully turbulent flow regime. 133 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 non- uniformity than fibers in the MD. 6.2.1 Further Comments 1) At high Reynolds numbers, the pressure drop correlates with the dynamic pressure , 2Vρ , and consequently the pressure coefficient with 0Re∝pC . In contrast, at low Reynolds number, the pressure drop correlates with viscous effects, LVμ , which leads to 1Re−∝pC . Therefore, considering the sum of the viscous and inertial resistance, the pressure drop may be expressed as 2VbVap ρμ +=Δ or 01 ReRe baC p += − . In this dissertation, pressure drop-flow rate relationships were of the form bp aC Re= . This agrees with some published work [2]. A more general form, however, is c p baC ReRe 1+= − [3]. 134 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. 135 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 sub- systems 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 μρ DU∞=Π1 (flow- cylinder property), dL=Π 2 and 33 LUEI ∞=Π μ (fiber and fiber-flow property), DL=Π 4 (fiber-cylinder geometrical property), cf=Π 5 (fiber-cylinder property), and Dχ=Π 6 , θ=Π 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. 136 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, =DL 2, 5, 10, and 20 oriented to the flow at angles between oo 900 ≤≤α and for 40Re1 ≤≤ 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 DL and Re , the drag coefficient increased monotonically with α , whereas the lift coefficient was approximately parabolic, achieving a maximum near o45=α . The curves of LC and DC against α were similar as one varies DL and Re . Both LC and DC decreased with increasing Re . 2) LC increased rapidly with DL for 10≤DL , and thereafter was nearly independent of DL . DC decreased rapidly with DL for 10≤DL , and thereafter was likewise nearly independent of DL . 137 3) Curve fits of ( )Re,, DLCL α and ( )Re,, DLCD α 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 o70≥α ), and overestimated the drag force by increasing amounts as α was diminished (error > 15% for o45=α and error > 50% for o30≤α ). 6) The non-dimensional torque coefficient is zero at inclination angles o0=α and o90=α . It should be noted that while o0=α is an unstable condition, a cylinder oriented with o90=α 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. 138 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. 139 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.