Applied Science, Faculty of
Mechanical Engineering, Department of
DSpace
UBCV
Shariati, Mohammad Reza
Doctor of Philosophy - PhD
University of British Columbia
The fluid flow patterns in a paper-machine headbox have a strong influence on the quality of the paper produced by the machine. Due to increasing demand for high quality paper there is a need to investigate the details of the fluid flow in the paper machine headbox. The objective of this thesis is to use experimental and computational methods of modeling the flow inside a typical headbox in order to evaluate and understand the mean flow patterns and turbulence created there. In particular, spatial variations of the mean flow and of the turbulence quantities and the turbulence generated secondary flows are studied. In addition to the flow inside the headbox, the flow leaving the slice is also modeled both experimentally and computationally. A Plexiglas scaled-down model of a headbox was constructed and LDV measurements of mean and turbulence velocities in the model were made. The numerical models used in Computational Fluid Dynamics (CFD) simulations were based on finite volume Reynolds-averaged equations and included turbulence representations with the standard k-e model, quadratic and cubic algebraic turbulence models and the Reynolds stress model (RSM). The "volume of fluid" method was used to model the free surface of the jet leaving the slice. The CFD investigation showed that the turbulence generated secondary flows in the headbox are much smaller than previously thought and that they are confined to a very small region close to the wall. Comparison of the experimental and numerical results indicated that streamwise mean components of the velocities in the headbox are predicted well by all the turbulence models considered in this study. However, the standard k-e model and the algebraic turbulence models fail to predict the turbulence quantities accurately. Standard k-e model also fails to predict the direction and magnitude of the secondary flows. Inaccuracies in the numerical prediction of turbulence quantities increase as the contraction ratio increases. This is due to the incorrect representation of the production of the turbulence kinetic energy by the k- e model and the algebraic turbulence models. Significant improvements in the k-e model predictions were achieved when the turbulence production term was artificially set to zero. This is justified by observations of the turbulent velocities from the experiments and by a consideration of the form of the kinetic energy equation. A better estimation of the Reynolds normal stress distribution and the degree of anisotropy of turbulence was achieved using the Reynolds stress turbulence model. Integral length scales were measured from the LDV velocity observations in the plexiglass scale model and were compared to the length scales obtained computationally using both standard k-£ and Reynolds stress turbulence models. Although similar magnitudes were found for the measured and computed length scales, different trends between numerical and experimental results were observed. Measured integral length scale and computed length scales have been reported to be similar only in decaying grid turbulence, so that the differences observed in this rapidly distorting flow are not surprising. The trends measured here, and their contrast to the calculated values of length scale, should be kept in mind when discussing turbulence characteristics in a rapidly converging section such as the headbox contraction. Careful examination of the measured turbulence velocity results shows that after the initial decay of the turbulence in the headbox, there is a short region close to the exit, but inside the headbox, where the turbulent kinetic energy actually increases as a result of the distortion imposed by the contraction. The turbulence energy quickly resumes its decay in the free jet after the headbox. The turbulence quantities obtained using the k-e model in the free jet have the same trend as the results from RSM model. However, the results using the k-e model are at unrealistically high compared to those found from the RSM turbulence model in the free jet. This is due to the error that is made in the k-e results from computations inside the converging section of the headbox. The RSM model gives more realistic values for the turbulence quantities in the free jet region when the R S M model is used throughout the contraction region as well. The overall conclusion from this thesis, obtained by comparison of experimental and computational simulations of the flow in a headbox, is that numerical simulations show great promise for predictions of headbox flows. Mean velocities and turbulence characteristics can now be predicted with fair accuracy by careful use of specialized turbulence models. Standard engineering turbulence models, such as the k - e model and its immediate relatives, should not be used to estimate the turbulence quantities essential for predicting pulp fiber dispersion within the contracting region and free jet of a headbox, particularly when the overall contraction ratio is greater than about five.
14554458 bytes
application/pdf
EXPERIMENTAL AND MATHEMATICAL MODELING OF FLOW IN HEADBOXES By Mohammad Reza Shariati B. Eng. D.M.E.T. College, Calcutta, India 1988. M.A.Sc. (Mechanical Engineering) Technical University of Nova Scotia 1995. A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIRMENT FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF MECHANICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA 2002 © Mohammad Reza Shariati In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mechanical Engineering The University of British Columbia 2324 Main Mall Vancouver, BC Canada V6T 1Z4 Date: March 2003 Abstract The fluid flow patterns in a paper-machine headbox have a strong influence on the quality of the paper produced by the machine. Due to increasing demand for high quality paper there is a need to investigate the details of the fluid flow in the paper machine headbox. The objective of this thesis is to use experimental and computational methods of modeling the flow inside a typical headbox in order to evaluate and understand the mean flow patterns and turbulence created there. In particular, spatial variations of the mean flow and of the turbulence quantities and the turbulence generated secondary flows are studied. In addition to the flow inside the headbox, the flow leaving the slice is also modeled both experimentally and computationally. A Plexiglas scaled-down model of a headbox was constructed and LDV measurements of mean and turbulence velocities in the model were made. The numerical models used in Computational Fluid Dynamics (CFD) simulations were based on finite volume Reynolds-averaged equations and included turbulence representations with the standard k-e model, quadratic and cubic algebraic turbulence models and the Reynolds stress model (RSM). The "volume of fluid" method was used to model the free surface of the jet leaving the slice. The CFD investigation showed that the turbulence generated secondary flows in the headbox are much smaller than previously thought and that they are confined to a very small region close to the wall. Comparison of the experimental and numerical results indicated that streamwise mean components of the velocities in the headbox are predicted well by all the turbulence models considered in this study. However, the standard k-e model and the algebraic turbulence models fail to predict the turbulence quantities accurately. Standard k-e model also fails to predict the direction and magnitude of the secondary flows. Inaccuracies in the numerical prediction of turbulence quantities increase as the contraction ratio increases. This is due to the incorrect representation of the production of the turbulence kinetic energy by the k- e model and the algebraic turbulence models. ii Significant improvements in the k-e model predictions were achieved when the turbulence production term was artificially set to zero. This is justified by observations of the turbulent velocities from the experiments and by a consideration of the form of the kinetic energy equation. A better estimation of the Reynolds normal stress distribution and the degree of anisotropy of turbulence was achieved using the Reynolds stress turbulence model. Integral length scales were measured from the L D V velocity observations in the plexiglass scale model and were compared to the length scales obtained computationally using both standard k-£ and Reynolds stress turbulence models. Although similar magnitudes were found for the measured and computed length scales, different trends between numerical and experimental results were observed. Measured integral length scale and computed length scales have been reported to be similar only in decaying grid turbulence, so that the differences observed in this rapidly distorting flow are not surprising. The trends measured here, and their contrast to the calculated values of length scale, should be kept in mind when discussing turbulence characteristics in a rapidly converging section such as the headbox contraction. Careful examination of the measured turbulence velocity results shows that after the initial decay of the turbulence in the headbox, there is a short region close to the exit, but inside the headbox, where the turbulent kinetic energy actually increases as a result of the distortion imposed by the contraction. The turbulence energy quickly resumes its decay in the free jet after the headbox. The turbulence quantities obtained using the k-e model in the free jet have the same trend as the results from R S M model. However, the results using the k-e model are at unrealistically high compared to those found from the R S M turbulence model in the free jet. This is due to the error that is made in the k-e results from computations inside the converging section of the headbox. The R S M model gives more realistic values for the turbulence quantities in the free jet region when the R S M model is used throughout the contraction region as well. The overall conclusion from this thesis, obtained by comparison of experimental and computational simulations of the flow in a headbox, is that numerical simulations show great promise for predictions of headbox flows. Mean velocities and turbulence characteristics can now be predicted with fair accuracy by careful use of specialized turbulence models. Standard engineering turbulence models, such as the k - e model and iii its immediate relatives, should not be used to estimate the turbulence quantities essential for predicting pulp fiber dispersion within the contracting region and free jet of a headbox, particularly when the overall contraction ratio is greater than about five. iv Table Contents Abstract ii Table Content v List of Tables vii List of Figures viii Acknowledgment xii 1 - Introduction 1 1.1 Headbox 4 1.2 Literature review 6 1.2.1 Distributor 6 1.2.2 Evening section and the slice 7 1.2.3 Free surface jet 9 1.3 Scope and objectives of this thesis work 10 2 -Computational Methods 12 2.1 Introduction 12 2.2 Governing Equations 12 2.3 Solution in Curvilinear Coordinates 13 2.4 Turbulence Modeling 14 2.5 Standard k-e model 14 2.5.1 Near-Wall treatment 15 2.6 Non linear Turbulence model 17 2.6.1 Quadratic Model 17 2.6.2 Cubic Model 18 2.7 Implementation of the Non linear Turbulence Models 19 2.8 The Reynolds Stress Model (RSM) 20 2.9 Free surface calculation using Volume of fluid methods 22 2.10 Turbulent Flow Calculation in Square Duct 23 3 -Experimental Arrangement and Measurement Techniques 29 3.1 Experimental Apparatus and Equipment 29 3.2 Laser Doppler Velocimetry 33 3.2.1 Background 33 3.2.2 Apparatus 34 3.2.3 Data Acquisition 36 v 3.2.4 Alignment 36 3.3 LDV operation parameters 37 3.4 Seeding particles 38 3.5 LDV Data Analysis 40 3.5.1 Removal of Bad Data Points 40 3.5.2 Calculation of Expected Values - Velocity Bias Correction 40 4 - Experimental measurements 43 4.1 Error analysis 43 4.1.1 Angular alignment of the LDV system 44 4.1.2 Measurement volume positions 44 4.1.3 Bulk velocity measurements 45 4.2 Velocity measurements results 47 5 - Numerical Simulations 62 5.1 Turbulent flow through the converging section of a headbox 62 5.1.1 Effects of the different boundary condition 63 5.1.2 Grid independence 76 5.1.3 Secondary flow calculations 79 5.2 -Flow through the headbox model 84 5.3 Turbulence Kinetic Energy 93 5.4 Turbulence Length Scale 101 5.5 Free surface Jet 105 6 - Summary and conclusion 112 6.1 Summary 112 6.1.1 Experimental measurements 112 6.1.2 Numerical simulation of the converging section 113 6.2 Conclusion 118 7 - Recommendation for future work 120 Nomenclature 121 Bibliography 124 Appendix A 130 Appendix B 133 vi List of Tables Table 2-1: Empirical constants used in the cubic nonlinear turbulence model 19 Table 2-2: Grid sizes for testing the grid independence of the square duct flow simulation 25 Table 3-1. The Geometry of the Headbox Converging Section 31 Table 3-2 : L D V beam properties 35 Table 3-3: Properties of the metallic coated hollow glass spheres 39 Table 4-1: Bulk velocity through different diffusers at the bottom row 46 Table 5-1: Computational geometry basic dimensions 62 Table 5-2: Boundary condition for numerical calculations 65 Table 5-3: Grid sizes for testing the grid independence of the flow simulation in the converging section of a headbox 76 Table 5-4: Boundary condition for numerical calculations 80 vii List of Figures Figure 1-1: Schematics of a Paper machine 3 Figure 1-2: Hydraulic Headbox 3 Figure 2-1: Vector plot of the secondary flows in the square duct 24 Figure 2-2Secondary flow calculation using the quadratic algebraic turbulence model for three different grid sizes at X=80x h along the Y-axis for Red=250000 (D h is the hydraulic diameter and Rej is the Reynolds number based on Dh). Y * =Y/a where a is the half width of the duct 25 Figure 2-3: Secondary flow calculation using the cubic algebraic turbulence model for three different grid sizes at X=80x h along the Y-axis for Red=250000 (D h is the hydraulic diameter and Rea is the Reynolds number based on Dh). Y * =Y/a where a is the half width of the duct 26 Figure 2-4 : Experimental and numerical calculation of secondary flows along the Y axis at X=80xDh and Red=250000 (D h is the hydraulic diameter and Re d is the Reynolds number based on Dh). Y * =Y/a 27 Figure 2-5: Secondary flow calculations along the Y axis at different distances from the entrance (quadratic nonlinear turbulence model). Y*=Y/a 28 Figure 3-1: The Flow loop in the experiment 31 Figure 3-2: Experimental model drawing. A l l dimensions are in meters 32 Figure 3-3 : Drawing of a typical diffuser. A l l dimensions are in meters. Wall thickness ~ 1 mm 33 Figure 3-4: L D V system schematics 35 Figure 4-1: Mean velocity profile measurement across different diffusers. Data are measured at the centerline of the diffusers on the bottom row. L * is the distance from the centerline of the center diffuser at the bottom row divided by the width of a diffuser 46 Figure 4-2 : Flow pattern inside a diffuser ( from White [51]) 48 Figure 4-3: Flow measurements across the width of the diffuser. Measurements are done 30 mm downstream of the exit plane of the diffuser. L * is the distance from the diffuser side wall nondimensionalized with the width of the diffuser. U is the streamwise component of the velocity 49 Figure 4-4: Streamwise velocity profile measurements at the centerline of a diffuser column. Measurements done 30 mm downstream of diffuser exit plane. Z* is the height from the bottom plane of the headbox divided by the total height of the headbox at that position 49 Figure 4-5: Measurement location along the center line of the symmetry plane 50 Figure 4-6: Mean velocity components measured experimentally in the region downstream the diffusers, along the line A B C in figure 4.5. Velocities are nondimensionalized with the bulk velocity at the inlet to the converging section. X * =X/Xi where X is the total length after the diffuser tubes 51 Figure 4-7: Rms. values of the fluctuating component of velocities u', v', w' normalized by U (local mean velocity in the X direction) along the line A B C in figure 4.5 52 viii Figure 4-8: Experimental measurements of rms. values of the fluctuating component of velocities u ' , v ' and w ' normalized by u 1 inlet along the line A B C in figure 4.5. 53 Figure 4-9: Mean component of the velocity in the streamwise direction nondimensionalized with the inlet bulk velocity. X * is the distance along the X axis normalized with the total length of the converging section 54 Figure 4-10: Rms. values of the fluctuating component of velocities u 1 , v 1 , normalized by u ' inlet 55 Figure 4-11 Schematic drawing of the velocity profile measurements location. Measurements are done along the vertical line at the center of a diffuser set (line AB) in the X direction 56 Figure 4-12: Mean velocity component in the X direction measured along the line A B in figure 4.9 at different X position. Z* is the distance from the bottom plane normalized with total height at the same location 57 Figure 4-13: Variation of the velocity profile in the Z direction at different location along the X axis. Z* is the distance from the bottom plane normalized with total height at the same location. Measurements are along the line A B in figure 4.9 58 Figure 4-14: Fluctuating component of velocity in the X direction (MD) measured along the line A B in figure 4.9 at different X position. Z* is the distance from the bottom plane normalized with total height at the same location 59 Figure 4-15: Fluctuating component of velocity in the Y direction(CD) measured along the line A B in figure 4.9 at different X position. Z* is the distance from the bottom plane normalized with total height at the same location 60 Figure 4-16: Fluctuating component of velocity in the Z direction measured along the line A B C in figure 4.5 at different X positions. Z* is the distance from the bottom plane normalized with total height at the same location 61 Figure 5-1: Computational geometry. A l l dimensions are in meters. 63 Figure 5-2: Geometry of the inlet for the case of the non-uniform inlet flow 66 Figure 5-3 : Computational mesh at the symmetry plane 66 Figure 5-4: Mean streamwise component of the velocity at the symmetry plane Y=0 for two different boundary condition at the exit with the uniform flow at the inlet along the line joining the center of the inlet plane to the center of the exit plane. U * is the mean streamwise component of the velocity divided by the bulk inlet velocity. X * is the distance in the X direction divided by the total length in the X direction 67 Figure 5-5: Mean component of the velocity in the Z direction at the symmetry plane Y=0. for two different boundary condition at the exit with the uniform flow at the inlet along the line joining the center of the inlet plane to the center of the exit plane. W* is the component of the velocity in the Z direction divided by the bulk inlet velocity. X * is the distance in the X direction divided by the total length in the X direction 67 Figure 5-6: Mean component of the velocity in the Y direction at the center line of the exit plane for two different boundary condition at the exit with the uniform flow at the inlet. V * is the mean component of the velocity in the Y direction divided by the bulk inlet velocity. Y * is the distance in the Y direction divided by the total length in the Y direction 68 ix Figure 5-7: Mean component of the velocity in the Y direction at the center line of the exit plane (Line A B above) for the case of the non-uniform flow at the inlet. V * is the mean component of the velocity in the Y direction divided by the bulk inlet velocity. Y * is the distance in the Y direction from the center of the line A B divided by the total length in the Y direction 69 Figure 5-8 : Cross flow contour plot for the non symmetric converging section with uniform inlet condition 71 Figure 5-9: Cross flow contour plot for the symmetric converging section with uniform inlet condition 72 Figure 5-10: Cross flow contour plot for the non-symmetric converging section with non uniform inlet condition 73 Figure 5-11: Cork screw motion in the non-symmetric convergence 74 Figure 5-12: Cross flow contour plot for the symmetric converging section with non-uniform inlet condition 75 Figure 5-13: Comparison of the mean component of velocity at the centerline for different grids. U*=U/U i n , X*=X/L X (see table (4.3 and 4.1) for definition of terms). Results obtained using R S M 77 Figure 5-14: Comparison of the turbulence kinetic energy, k* = k/kjn. X*=X/Lx , Here kj n is the turbulent kinetic energy at the inlet to the headbox. Results obtained using R S M 78 Figure 5-15: Comparison of the CD velocity at the exit plane, V*=V/Ui„ , Y * = Y / L Z , U i n is the bulk velocity in the streamwise direction at the inlet. Results obtained using R S M 78 Figure 5-16: Mean component of velocity at the centerline. U*=U/Uj n, X*=X/L X (see table (4.3 and 4.1) for definition of terms) 81 Figure 5-17: CD velocity at the exit plane, V*=V/U i n , Y * = Y / L Z , U i n is the bulk velocity in the streamwise direction at the inlet 82 Figure 5-18: Turbulence Kinetic energy, k* = k/k i n. X*=X/Lx , Here kj„ is the turbulent kinetic energy at the inlet to the headbox 83 Figure 5-19 : Computational geometry 85 Figure 5-20: Computational grid 86 Figure 5-21: Computational grid details in the wake region of the diffuser walls 86 Figure 5-22: Comparison of the mean velocity components in X direction in the region after the diffusers for different grid sizes. Velocities are nondimensionalized with the bulk velocity at the inlet to the converging section. X * =X/Xi where X i is the total length after the diffuser tubes 87 Figure 5-23: Turbulence kinetic energy variations in the region after the diffusers for different grid sizes, k* is the turbulence kinetic energy nondimensionalized by the 3/2 (Bulk velocity)2 at the inlet to the converging section 88 Figure 5-24: Streamwise component of the velocity measured along the vertical line at the center of the diffusers at the location 3 cm after the diffuser exit plane. Z* is the height from the bottom plane divided by the total height at the location 90 Figure 5-25: RMS values of the fluctuating component of the velocity in the streamwise direction measured along the vertical line at the center of the diffusers at the location 3 cm after the diffuser exit plane. Z* is the height from the bottom plane divided by the total height at the location 91 Figure 5-26: Comparison of the mean velocity components in the region after the diffusers. Velocities are nondimensionalized with the bulk velocity at the inlet to the converging section. X * =XJX\ where Xi is the total length after the diffuser tubes . 92 Figure 5-27: Turbulence kinetic energy variations in the region after the diffusers. k* is the turbulence kinetic energy nondimensionalized by the measured k at the inlet to the converging section located at x =0.6185m(6 mm after the convergence starts).. 93 Figure 5-28 : Rms. values of the fluctuating component of velocities u 1 , v 1 and w ' normalized by u 1 inlet 95 Figure 5-29: Different components of the Reynolds shear stresses normalized by u'u' at inlet using the Reynolds stress model 96 Figure 5-30: Turbulence Kinetic energy at mid plane converging section 97 Figure 5-31: Comparison of the Production term and rate of dissipation term(epsilon) in the k equation for different models CR=10 97 Figure 5-32: Geometry of a 2 dimensional converging section 99 Figure 5-33: Turbulence kinetic energy k calculated for different contraction ratios.... 101 Figure 5-34: Variation of the length scale along the centerline in the mid plane near the entrance. L * is the length scale divided by the width(=0.0375m) of the diffusers at the exit. X * is the distance in the x direction after the diffuser tubes divided by the total length of the headbox after the diffuser tubes 104 Figure 5-35: Free surface jet calculation domain 108 Figure 5-36: Numerical Grid used for free surface calculation 109 Figure 5-37: Mean component of the velocity in the streamwise direction nondimensionalized with the inlet bulk velocity. X * is the distance along the X axis normalized with the total length of the converging section 109 Figure 5-38: Rms. values of the fluctuating component of velocities u ' , v 1 ,w 1 normalized by u 1 inlet 110 Figure 5-39: Turbulence kinetic energy along the centerline of the mid plane inside the headbox and along the center of the free surface jet outside the headbox. X * =1 is the exit -plane of the headbox. K * is the turbulence kinetic energy normalized with k value at the inlet I l l Figure 5-40: Variation of Reynolds shear stresses in the free surface jet. Reynolds shear stresses are normalized with 2/3 k at the inlet. X * is the distance along the X axis divided by the total length of the headbox I l l xi Acknowledgments I would like to express my deepest gratitude to my research supervisors Dr. Martha Salcudean and Dr. Ian Gartshore for their continuous guidance, inspiration and their endless support. Appreciation is due to Dr. R. J. Kerekes, Dr. Calisal and Dr. Dan Fraser for participating on my research committee and for their help and guidance. I would like to thank all the people who helped me in my research. I owe a great deal to my colleagues at U B C specially Mike Georgallies, Suqin Dong, and Dr. Eric Bibeau. Special thank to Dr. A l i Mahallati, Mr. Kourosh Baseghi, Mr. Fari Fouladi, Mr. Babak Samei and Mr. Bardia Goodarzi for their endless help and support in my research work and for being a true friend since I have known them. I also wish to acknowledge the financial assistant from Forest Research of B C and Natural Sciences and Engineering Research Council. Finally, I am forever indebted to my wife, Jennifer, my kids, Emily, Hannah, Farah and Keyan, my parents and my brothers and sister for their understanding, endless patience and encouragement when it was most required. xii 1 - Introduction The art of true paper making had its origin in China around 100A.D. It was discovered that filtering a pulp suspension through a fine screen could form a mat of fiber. The mat after pressing and drying was found to be suitable for writing and drawing. Prior to the industrial revolution, papermaking produced a useful product but suffered from size and output limitations. However, at the end of the 18th century the process of papermaking made a major step forward in terms of output with the introduction of continuous forming. The Fourdrinier brothers in Britain introduced a practical continuous former in 1810 following development work carried out in France. The early machines consisted of a flowbox adding paper stock to the moving wires supported between two rolls. While the basic steps in the process are essentially the same as those used in the early craft, papermaking has made many advances since the introduction of continuous forming and has become an advanced industry. Modern papermaking can be divided into two main parts: the first part is the processing of the pulp and the second part is the formation of the actual paper in the paper machine. After being processed, a blend of one of more types of pulp with fillers and additives (generally called stock), is pumped to the paper machine's headbox. The function of the headbox is to distribute the stock evenly and to control alignment of the fibers in both the machine direction and across the width of the paper web (the cross direction) as the stock is fed onto the forming section through a slice opening. Figure 1-1 and 1-2 show the role of the headbox in papermaking and the function of a typical hydraulic headbox, respectively. The forming section performs the same function as the screen used by early papermakers. This part of the machine drains water from the stock, leaving behind the fibers which naturally bond together into a sheet. The Fourdrinier forming section features a continuous wire loop that allows the water to drain from the stock. Other more recent designs include the twin wire former, which eliminates the two-sidedness problem associated with the Fourdrinier (i.e. visual differences between the two sides of the paper sheet). The twin wire system l simultaneously forms the top and bottom of the web as the fibers are squeezed between the wires. The press section follows the forming section of the paper machine. In this part of the machine, 'nips' (large rollers) squeeze the water from the web. The web is supported on a looped felt mat (the paper machine felt) as it passes through the press section. At this stage after leaving the press section, the web will still contains up to 60% water (which means the same as a consistency of 40%). The area of the paper machine before the press section is called the wet end, while the remainder of the line from the press section onwards is called the dry end. After the press section, the web proceeds to the drying section, which consists of a series of steam heated cylinders configured into several groups, and often enclosed in a hood structure. The paper web dries by evaporation as its passes over the cylinders, or between the cylinders in the case of two tier dryers. Following drying, the paper web then passes through a calender. The calender is comprised of heated steel rolls that help to give the web smoothness and increase its density. The line is then completed by the winder, which winds the paper web into jumbo rolls that can be two meters or more in diameter. Rather than being a single unit, a paper machine is in reality a complicated series of individual machines designed to operate as a total unit, with a sheet of paper as a common thread. A critical element in this complex chain is the headbox. Each section, or unit, of the paper machine is highly dependent on the unit preceding it. The headbox and flow approach system unit is the most important section of this multi-machine complex. If the flow from the headbox is not correct then the Fourdrinier and subsequent machine sections cannot function properly [1]. In the past, due to a lack of suitable models for the flow in the headbox, most improvements in headbox design were made by trial and error. In recent years, the increase in computational power and advances in the CFD (computational fluid dynamics) field have enabled researchers to model the headbox flow numerically. However, the complexity of the flow coupled with its enormous size have forced researchers to make simplifying assumptions in modeling the headboxes numerically. 2 These simplifications and a lack of reliable experimental data have led to some questions regarding reliability and usefulness of these simplified numerical solutions. To answer some of the questions the present numerical and experimental modeling of the headbox flow was undertaken. - O x . Presses H e a d B o x Evaporation u Drainage Paper Dryers Dewatering Figure 1-1: Schematics of a Paper machine Figure 1-2: Hydraulic Headbox 3 1.1 Headbox The function of the headbox is to take the stock delivered by the pump and transform the pipe flow into a uniform rectangular jet equal in width to the paper machine and at uniform velocity in the machine direction. Since the formation and uniformity of the final paper are dependent on the even dispersion of the fibers and fillers, the design and operation of the headbox system is absolutely critical to a successful paper making system. The following may be listed as primary objectives of the headbox system: 1. To accelerate the stock to a speed close to that of the wire formers. 2. To spread stock evenly across the width of the machine 3. To remove or minimize cross-machine consistency variations. 4. To remove or minimize machine direction velocity gradients. 5. To create a level of turbulence which will eliminate fiber flocculation. 6. To discharge an even flow from the slice opening so that the stock will impinge on the wire at the correct location and angle. 7. To control the fiber orientation in the stock leaving the slice so that fibers are distributed uniformly and with a known (preferably controllable) orientation to create the required properties in the resulting paper. In addition to these requirements, there are some secondary ones as well as some economical ones. A detailed discussions of an ideal headbox requirements is given in [2] There are three basic types of headbox design [3]: 1. Open headboxes 2. Air-cushioned pressurized headboxes. 3. Hydraulic headboxes Open headboxes were employed in the early paper machines, where the gravity head of stock was used to give the correct discharge velocity. As the speed of paper machines increased, it became impractical to further increase the height of the stock, and the pressurized headbox was developed. These headboxes depend on the action of the rotating perforated rolls to provide turbulence and level out velocity gradients. The main drawbacks of air-cushioned headboxes are: • Air bubble entrapment in the stock. 4 • Use of the moving parts to generate turbulence for dispersion and deflocculation of the fibers. The high speed and special characteristics of modern paper machines demanded better headboxes and this has led to the development and wide spread use of hydraulic headboxes. A review of different types of headboxes and their applications is given in [3]. Elimination of the free surface in hydraulic headboxes allowed the use of stationary elements for flow evening and turbulence generation at much higher flow velocities. The higher velocities and use of stationary elements have radically increased turbulence intensity levels, which has improved deflocculation. Because of their improved performance and competitive price, hydraulic headboxes are now the most common type in operation in today's paper mills [4]. All hydraulic headboxes consist of three major parts [5]: the distributor, the evening section and the slice. The distributor deals with the initial distribution of the flow from the pipe to the full width of the machine. The flow from the distributor is improved through the rectification and evening processes. In a hydraulic headbox, a tube bank is often used in these processes. The wall friction in the tubes dampens flow non-uniformities originating in the stock approach system and creates turbulence, which is needed to prevent fiber flocculation in the paper-machine forming zone. The processes in the tube bank may include mixing and blending of separate flows from a distributor, eliminating undesirable cross-flow and eddies, improving the velocity profile and developing turbulence of desired scale and intensity. The jet development process is done in the slice and can be described as delivering the stock to the sheet forming section at a speed close to the wire speed. An ideal headbox should produce a uniform and stable jet over the width of the machine, without lateral velocities and machine-direction perturbations. The headbox flow, particularly the turbulence levels in the headbox, has a profound effect on the fiber orientation and ultimately on the finished paper quality [5]. Therefore it is necessary for the headbox to generate and maintain a certain degree of isotropy and scale of turbulence. Although all these three parts are integrated and considered as one unit, past researchers have often studied the flow in the different parts of the headbox separately. 5 1.2 Literature review 1.2.1 Distributor How well the headbox can spread the stock across the width of the paper machine is dependant on the design of the flow approach system. Early designs of headboxes accepted feed flows from either single or multiple entry points, combined them where necessary, spread them across the width of the machine, and directed the flow towards the slice. Such designs frequently fell short of providing adequate flow uniformity, particularly as machine speeds increased and the demands of cross-machine paper uniformity became more stringent. Superior results have been achieved in modern headboxes by employing a manifold distributor feeding into an array of relatively small diameter tubes. The advantages of this system are that it is simple in design, compact, relatively easy to construct and flexible in operation [6]. The analysis and design of headbox manifolds was traditionally based on a one-dimensional energy balance of the flow inside the manifold [7,8]. A detailed consideration of the principles of the tapered manifold type flowspreaders is given in [9] where design guidelines to obtain a uniform flow distribution at the slice exit can be found. The early numerical modeling of the headbox manifold was done with very simplified geometrical models of the headbox and could not show flow details. Later numerical modeling was introduced as an additional tool in the study of headboxes. Among such investigations is the work done by Syrjala et al. [10] who modeled an actual headbox manifold without the diffuser tubes. They compared their results with measurements done in air and concluded that the numerical results can be used reliably. Jones and Ginnow [11] used a commercial computational fluid dynamic (CFD) software (Fluent) employing the k-e turbulence model to predict turbulent flow in several types of flow regime encountered in the headbox in simplified form. They calculated the flow distribution in a straight channel diffuser in three dimensions and in a Beloit experimental headbox in two dimensions. Predictions compared favorably with available experimental data but required an arbitrary change in some of the constants used in the turbulence model. Their main conclusion was that reliable experimental data were needed for 6 verification and tuning of the turbulence models. The authors recommended further experimental investigation of the eddy size variation in manifolds. Shimizu and Wada [12] calculated a generic headbox using a body-fitted coordinate system. The flow distribution in the manifold was investigated in two dimensions and the flow area reduction was calculated in three dimensions with assumptions of periodicity. The jets from the diffuser tubes were modeled in three dimensions using calculation results obtained for a single tube. Hamalainen [13] linked two-dimensional models of the manifold, the turbulence generating section, and the slice using a finite element method. The pressure drop in each tube was assumed to be that of a single tube using the homogenization technique. The model was then applied to the optimization of a headbox [14]. Separate three-dimensional models of the manifold and of the slice have been performed by Lee and Pantaleo [15] to investigate various effects of headbox control devices on flow characteristics and fiber orientation. They verified the important effect of the recirculation on the inlet velocity profile. Most of these works studied the manifold and flowspreader in a simplified form and in isolation from the rest of the headbox, the main simplification being the replacement of the tube bank with a slot or, in a few cases, replacing it with the results of a flow calculation for a single diffuser tube. This results in a very idealized flow, different from that in an actual headbox. A complete modeling of the manifold along with all the individual tubes was done by He et. al.[16] and Hua et. al.[6,17,18,19]. They concluded that good distribution of flow could be maintained over a wide flow range by using appropriate values of recirculation. 1.2.2 Evening section and the slice Like the manifold section, the analysis and design of the evening section and slice were traditionally based on simple one-dimensional flow models. Until recently, design improvements have been based on trial and error. Among the early studies are the work done by Lee and Pantaleo [15]. They studied the effects of the slice lip profiling on the ratio of the cross machine velocity (CD) to the machine direction velocity (MD). Among the more recent works is the study of the Aidun and Kovacs [20] in which the converging 7 section of a commercial headbox was modeled. These authors showed a large cross flow at the exit from the headbox. The classic generation of secondary flows in a square duct due to anisotropy of the turbulence, led them to believe that significant secondary flow is generated in the converging section of the headbox. They also believed that these secondary flows have a large magnitude and extend over a large area in the exit plane, therefore having a significant adverse effect on the fiber orientation and paper quality. They also concluded that simple turbulence models such as standard k-e are inadequate for modeling these secondary flows. In the follow-up study of Aidun [21] the existence of these secondary flows was shown, and their effect on the fiber orientation was also discussed. Different turbulence models namely R N G k-e and Reynolds stress model for modeling the headbox flow for different contraction ratios (CR) were compared by Bandhakavi and Aidun [22]. The authors found the results from the Reynolds stress model superior to that of the R N G k-e model. However, they suggested further study of the turbulence convective transport and the stress production term. Parsheh and Dahlkild [23] modeled the mixing of the different layers of a three-layer stratified headbox jet using commercial fluid dynamic software ( C F X ). Their main conclusion was that vane length has a considerable effect on mixing while the effects of the shape of the vane tip has a minimal effect on mixing. The performance of different turbulence models for computing the turbulent flow at the centerline of a two dimensional contraction was investigated by Parshe [24]. The results were compared to measured data from a two-dimensional contraction which was tested in a wind tunnel. It was observed that the numerical prediction of the turbulence development in the converging section of the headbox by different models was dependent on the prescribed inlet boundary condition. The authors also reported that the turbulence kinetic energy computed by the k-e model was exaggerated. Previous numerical studies for the calculation of the flow distribution from the manifold either ignore the diffuser tubes, model tubes in two-dimensions, or assume single tube behavior. These studies may not model the flow non-uniformities existing across the diffuser tubes caused by 1. Headbox design problems and design compromises resulting in deviation from the desired pressure profile across the manifold. 8 2. Variation in geometry of the tubes 3. Change in the pressure drop across identical tubes caused by velocity variations entering the tubes produced by the manifold. Numerical modeling of the slice section does not in general incorporate the exact flow field emerging from the turbulence generating section. A detailed three-dimensional numerical modeling of the manifold along with the modeling of each individual tube in the tube bank and the converging section of an industrial size headbox was performed by Hua et al [17,18,19] using the standard k-e turbulence model. The major finding of their study was that the flow at the inlet of the slice contains local flow non-uniformities and flow variations in the paper thickness direction (PD) and the cross-machine direction, creating flow structures at the beginning of the slice. They also showed that mean flow structures at the inlet caused by the footprint of individual tubes would disappear towards the exit. However any cross machine nonuniformity originating from a manifold flow imbalance, would still exist at the exit. The headbox flow effect on fiber orientation was studied by Zhang et. al [5,25] experimentally and numerically. They measured the orientation of dyed nylon fibers moving in a pilot plexiglass headbox and compared with the numerical results obtained by simulation of the fiber motions in the headbox. Zhang concluded that the anisotropy of the fiber orientation in the headbox flow is caused not only by the mean flow field characteristics, but also by the turbulence characteristics. Therefore for accurate modeling of the fiber motion in the headbox the explicit effects of the turbulence must be considered. 1.2.3 Free surface jet The free surface jet emanating from the headbox has a profound effect on the fiber orientation, and therefore on the final paper quality. Few headbox simulations include the free surface jet, but many experimental studies have been done in the past. Studies of the jet from the headbox can be divided roughly into three categories: fiber motion in the jet, surface streaks and instability and effects of the jet angle and speed on the dewatering at the impinging zone. 9 Ullmar and Norman [26,27] made experimental observations of fiber orientation in the headbox jet. In their studies a CCD- camera was mounted over the outlet of a symmetric headbox model to observe the fibers. The experiments were done for different contraction ratios and different flow rates. The result of their study showed a strong relation between contraction ratio and fiber orientation. They also indicated a very small effect of the flow rate on the fiber orientation. A high speed digital imaging of paper forming technique was developed by Aidun and Ferrier [28] and was used to investigate the origins of the streaks in the jet and in the early forming zone [29]. They speculated that the interaction of jets from two rows of tubes generates a mean secondary flow of counter-rotating vortices for each pair of jets exiting from the diffuser tubes and these vortices generate streaks at the slice and on the forming table. However, Soderberg [30,31] concluded that the origin of the streaks was the disturbances inside the headbox. Soderberg also concluded that the most dominant streaky structure was connected to the break-up of spanwise waves. Soderberg believed that spanwise waves are the result of the velocity profile relaxation process and can be explained by hydrodynamic stability theory. In general, it appears that there is no consensus on the origin and cause of the streaks that appear in the jet and on the forming tables and further studies are needed. 1.3 Scope and objectives of this thesis work The current research is part of an effort at the University of British Columbia to develop a computational method as a general design tool for accurate simulation of the fluid flow and fiber motion in complete headbox. Such a computational tool could be used to correct or improve existing headboxes, evaluate proposed retrofits, compare headbox designs, predict the influence of control devices and operation modes on flow behavior. It can also be used to serve as an important step for predicting fluid-fiber interaction in a way that can be directly applicable and beneficial to the pulp and paper industry. The main objectives of this thesis can be listed as follows: 1. To study the flow in the converging section of a headbox both experimentally and numerically. 10 2. To provide an experimental data set, representative of the flow in the strongly converging section of the headbox, which can be used to understand the turbulence in this section. 3. To assess through comparison between experimental and computational data the capability of the numerical models to predict the flow. To accomplish these objectives a scaled plexiglass headbox was built at the Pulp and paper center located at the University of British Columbia, and used in the present experiments. Laser Doppler Velocimetry (LDV) was used to measure the mean and turbulence quantities in the flow at several locations in the model. Two new numerical turbulence models were added to the existing CFD code (CMGFD) developed by Dr. Paul Novak [32] at the University of British Columbia. The two new turbulence models were used to model the turbulence generated secondary flows in the headbox and the numerical results were compared with experimental values and the values obtained from FLUENT software (Version 5.3) using the Reynolds Stress Model (RSM) for the same headbox flow. Finally FLUENT software was used to study the flow in the headbox and the free surface jet after the headbox. 11 2 -Computational Methods 2.1 Introduction One of the objectives of this thesis is to develop a numerical method for the accurate prediction of the flow field in a headbox. The appropriate model should be able to predict the mean velocity field and turbulence characteristics with an accuracy, which reflects allowable physically plausible conditions. In this chapter, the governing equations used in the development of the numerical model are described, and several different turbulence models are discussed. 2.2 Governing Equations We consider here the for the steady state conditions turbulent flow of a viscous, incompressible fluid with constant properties, and no heat transfer. The governing field equations are the Navier-Stokes and continuity equations, which are given for constant density and viscosity by ^ at oxj J axj dXj I ^ dXj dx,. J d U . }-=0 (2) d x . J where Ui is the velocity vector, p is the pressure F is the body force(which can include a gravitational potential), and \i is the viscosity of the fluid. A usual time averaging process on equation 1 and 2 for the steady state condition, leads to the time averaged continuity equation and the Reynolds-averaged Navier-Stokes equation which are given, in Cartesian tensor notation as 12 d U j U i d x. (VSy-u'tUj) (3) = 0 (4) where f d U i U j + — — -d Xj (7 X; V and p w,*/,. = Tfj is the Reynolds-stress tensor. In these equations, the over-bar and prime indicate the time-averaged mean and the fluctuation from the mean value, respectively. 2.3 Solution in Curvilinear Coordinates In order to obtain a numerical solution to the governing equations presented in the previous sections, the equations are discretized using the finite volume method. The resulting system of algebraic equations is then solved iteratively. In the following sections, the solution method in curvilinear coordinates is presented. More details on the computational approach used here for general, curvilinear coordinate systems may be found in He and Salcudean [33], and He [34]. These methods have been implemented in a curvilinear, multigrid, finite volume code called C M G F D developed by Paul Novak in the Mechanical Engineering Department of the University of British Columbia [32]. In the C M G F D code, the physical geometric quantities for each control volume are used to formulate the numerical scheme. The transport of a general quantity (j) is broken into orthogonal and non-orthogonal components. The orthogonal component of the flux has the same form as for a Cartesian co-ordinate system and methods such as the hybrid scheme, power-law scheme or the 13 second order schemes for regular geometries are applied to these terms. A second order accurate scheme is used for calculation of the non-orthogonal terms. The physical tangential velocity components are used by the C M G F D code as the dependent variables for the solution in curvilinear coordinates. Also a staggered grid arrangement is adopted in which the pressure is located at the geometric center of the control volume and the tangential velocity components lie at the midpoints of the respective control volume surfaces [34]. 2.4 Turbulence Modeling The objective in turbulence modeling is to represent the unknown Reynolds stresses in terms of known parameters. One of two main categories of modeling approaches is used for most engineering problem solutions. The first, known as turbulent-viscosity or eddy viscosity modeling, is based on Boussinesq's suggestion [35] that the Reynolds stresses are proportional to the mean strain rates. The second category is based on the development of differential equations describing the transport of individual stresses. 2.5 Standard k-e model As a result of the simplicity of their structure, linear turbulence models based on the eddy viscosity hypothesis are most commonly used in the solution of practical engineering problems. Among linear turbulence models, the k-e model is the most widely applied. Since the linear k-e model has been calibrated with thin layer shear flows it usually does quite well in the description of non-separated turbulent boundary layers which are encountered in many engineering problems. The long historical development and success of the usual k-e model, along with the fact that it can be incorporated into most Navier-Stokes computer codes in a relatively simple manner, constitute the main reasons for the continued popularity of this approach. In the linear k-e model, the Reynolds stress tensor is given by _ 2 (5) Tij-vijSij~~^kSij 14 where k=l/2 u:u: and S!t = '' d Ui • + d U j A d x, d x ; (Repeated suffices imply summation over all directions). k is the turbulent kinetic energy, and v t denotes turbulent kinematic viscosity which, unlike its laminar counterpart, varies spatially, and is not a property of the fluid. This model introduces transport equations for the two turbulence quantities, namely, the turbulence kinetic energy k and its rate of dissipation £. The k equation is: 9 ( * U , ) = 9 ' <?X; dx v + -j V dk •kj dx, + P-£ (6) where P = -uiuJSu and£=vS^ P and e are the rates of kinetic energy production and dissipation per unit mass, respectively. The transport equation for the kinetic energy dissipation is given in the same form and is: d -e U dx. J J dx. J V + - t d 6 6 +(C P - C e ) -dx. 1 2 k J (7) The isotropic viscosity assumption for the kinematic viscosity in terms of k and 6, given by v t= k 2 /e, finally closes the system of equations 3,4,6, and 7. The usual values of the constants are[36,37] C t x = 0.09, d = 1.44 , C2=1.92 , Gk=1.0 o e = K 2 / [ ( C 2 - C , ) C / 2 ] Where K=0.41 is the Von Karman constant. 2.5.1 Near-Wall treatment The standard k-e model described in the section is valid only for the fully turbulent regions and does not take into account the viscous effects which are important in the 15 vicinity of solid boundaries. Close to the wall the turbulent Reynolds number becomes small and the viscous effects dominate over the turbulence effects. To take this effect into account without using an excessive number of grid points near the wall, the standard wall function suggested by Launder and Spalding [38] is adopted in the present study. The standard wall function is based on the assumption that the flow in the vicinity of the wall behaves locally similar to an equilibrium turbulent boundary layer flow. Therefore the velocity profile near the wall can be calculated based on the value of the dimensionless normal distance from the wall y + as follows: U n 1 —P- = Mn(Ey+) , i f / > 11.36, uT K —?- = y+ , i f / < 11.36 «T where the dimensionless distance / is defined as: + WpUT y =—-— and the friction velocity is defined as: where y p is the normal distance of the first grid point (p) from the wall, U p is the velocity component parallel to the wall at the node p, and k and E are empirical constants which are taken as K= 0.41 , E - 9.97 (for the smooth wall) following Launder and Spalding [37]. By knowing the U p and y p the friction velocity ux can be found by solving the above equations iteratively. The values of the turbulence kinetic energy k and its rate of dissipation e at the first node (p) from the wall are then specified by the following expressions: 16 The use of a standard wall function provides the boundary conditions at the solid wall for the turbulence kinetic energy, and its rate of dissipation and the wall shear stress x w . Detailed descriptions may be found in Launder and Spalding[36,37] and He and Salcudean[33]. 2.6 Non linear Turbulence model Eddy-viscosity models have a major problem associated with them: • They are oblivious to the presence of the rotational strains (e.g. they fail to distinguish between the physically distinct cases of plane shear, plane strain, and rotating plane shear). In order to decrease the effects of this deficiency, one has to use a nonlinear or anisotropic generalization of eddy-viscosity models. These models express the Reynolds-stress tensor as a higher order function (normally 2 n d or 3 r d order function) of the mean velocity gradient. 2.6.1 Quadratic Model One of the non linear k-e models used in this study is the one suggested by Speziale[39,40]. In this model, the Reynolds stress tensor is expressed as a quadratic function of mean velocity gradients. One of the advantages of this model is its similarity to a standard k-e model. It basically constitutes a two-equation model with an anisotropic eddy viscosity. The two equations are, as before, the transport equation of the turbulence kinetic energy k (Eq. 6) and dissipation rate e(Eq. 7). The shear stress can be calculated from 17 c _r r2 — 3 / . S;K S t; S m„ S S ik kj *y mn mn ij _ 2 C C 2 — L e * e2 f. 1 (8) where ~ _ dSy d Uj _ «?* +Ukaxk~dxks»~dxk Ski Values of the two constants are set to CD=C e =1.6 as suggested by Speziale[39]. 2.6.2 Cubic Model In the present study to assess the performance of the different models a cubic model was also implemented. The cubic model used in this study is that suggested by Launder[41]. In this model, shear stress can be calculated from; 9 v v T..=-kS..—t-s.. + a - L v 3 u k v 1 £ ^ik^kj 3 ^ kl^kl^ij + C, V {aikski+Qikski \ llx txl Jtx txl J £ [aikQik~\aklQklSii) \ llx J fx txl fxl y J + C..v f t p . p SJ.QJ.+SJ.Q,.—S, Q., 8.. ki Ij kj li 3 km Im ij Skl\ + C6SijSklSkl+C7SijQkPkl 18 s..= y rdU. dU? I , J dx. J + dx. I) Q... = V i J dX. d X.\ (12) where C\ through C7 are empirical constants, and their values as suggested by Launder[41] are given in Table 2-1. Table 2-1: Empirical constants used in the cubic nonlinear turbulence model. Ci c 2 c 3 c 4 c 5 c 6 c 7 -0.16 0.1 0.26 -10CV2 0 -5C^ 2 5C^ 2 A drawback of this model compared to the quadratic model is that it requires a larger number of empirical constants. 2.7 Implementation of the Non linear Turbulence Models The C M G F D code with the standard k-e uses a linear relation between Tjj and the rate of strain tensor Sy. This approximation results in having the same form as the viscous term in the momentum equation. This similarity greatly simplifies the discretization of the momentum equations. To add the capability of solving turbulent flows using the nonlinear k-e models to the existing code, one has to calculate the Tjj for each momentum equation separately. Furthermore Tjj no longer has a similar form to the viscous terms due to the addition of the quadratic and cubic terms. Each of the terms in equations 8 and 12 has to be discretized separately and the contribution of each term to the diagonal and off diagonal terms in the solution matrix has to be added to the appropriate term in the solution matrix to obtain a final solution form. However, this requires changing the whole solution procedure for the momentum equation in the existing code. To avoid these complications the following simplifying assumptions are made: 19 1. Equations 8 and 12 are divided in two parts, a linear part and higher order terms. Since the linear part is already implemented in the code only the addition of the higher order terms is necessary. 2. Higher order terms are calculated separately and added to the equation as a source term. This greatly reduces the cumbersome task of discretizing each term and adding them to appropriate terms in the solution matrix. However, this may reduce the diagonal dominance of the solution matrix and consequently slows the convergence rate. 2.8 The Reynolds Stress Model (RSM) As pointed out by Speziale [39], the two-equation turbulence models based on the eddy viscosity assumption have the following major deficiencies: 1. The inability to properly account for stream line curvature, rotational strains, and other body force effects, aside from gravity. 2. The neglect of non-local and history effects on the Reynolds stress anisotropy. Some improvements can be achieved by the use of two equation models with nonlinear algebraic models, but to overcome the second of these deficiencies a higher order closure model such as a Reynolds stress model must be used. The Reynolds stress model involves calculation of the individual Reynolds stresses using differential transport equations. The individual Reynolds stresses are then used to obtain closure of the Reynolds-averaged momentum equation. In addition to the equations for the individual Reynolds stresses an extra equation is also required for turbulence dissipation, e. This brings the total number of additional equations to be solved for 3 dimensional computations to 7. By multiplying the instantaneous momentum equations with the fluctuating velocity, then taking the Reynolds-average of the product, one can derive the exact form of the Reynolds stress transport equations. Unfortunately, several of the terms in the exact equation are unknown and need to be either evaluated using the higher order moment equations (which themselves include further unknowns), or modeled in terms of other known variables. A typical equation for the Reynolds stress as described in FLUENT 5 user's guide, can be represented as 20 s V ' * - v ' Local Time Derivative Qj = Convection - ^ [p wujUk+p{SkjUi + 6ikUj)\ + ^ Dlf = Turbulent Diffusion = Molecular Diffusion V V ( DU' dU" \ - p \UiUk ^ + UjUk - pPigWjO + 9jut0\ P i j = Stress' Production ^ s P r o d u c t i o n <f>ij = Pressure Strain c*i = Dissipation Fij = Production by System Rotation Here £yk is the rotational third order tensor defined through the cross product of two vector as : Aj x Bj = Bp Ck. Of the various terms in these exact equations, Cy, DyL, Py, and Fy do not require any modeling. However, D/, Gij,^ , and e,y need to be modeled to close the equations. A commercial CFD Package FLUENT (version 5.3) was used to carry out the Reynolds stress model calculation in the present study. The Reynolds stress model employed by the FLUENT package is a second moment closure model and hence models the terms Ay , Gij,tyij , and ty in terms of known variables. More details for the modeling of this terms can be found in FLUENT 5 user's guide. 21 2.9 Free surface calculation using Volume of fluid methods The VOF model is a fixed grid technique designed for two or more immiscible fluids where the position of the interface between the fluids is of interest. In the VOF model, the volume fraction of each of the fluids in each computational cell is tracked throughout the domain. In each control volume, the volume fraction of all phases sums to unity. The properties appearing in the transport equations are determined by the presence of the component phases in each control volume. In a two-phase system, for example, i f the phases are represented by the subscripts 1 and 2, and i f the volume fraction of the second of these is being tracked, the density in each cell is given by p=a 2p 2 + ( l-a 2 ) pi where a,, is the volume fraction of the n t h fluid in the cell. A l l other properties (e.g., viscosity) are computed in this manner. The solution of a continuity equation for the volume fraction of one (for a two-phase system) of the phases determines the interface separating the two fluids. For the n t h phase this equation takes the following form: d(X„ dcXy, —n-+ui—^- = 0 dt dxj where i= 1,2,3. A single momentum equation similar to equation 1 is solved throughout the domain, and the resulting velocity field is shared among the phases. The momentum equation, is dependent on the volume fractions of all phases through the properties p and \i. In the case of the turbulence quantities, a single set of transport equations is solved, and the variables (k,e or the Reynolds stresses) are shared by the phases throughout the field. The VOF model in the CFD Package FLUENT was used for tracking the free surface jet which exits from the headbox. More details for the VOF modeling can be found in F L U E N T 5 user's Guide. 22 2.10 Turbulent Flow Calculation in Square Duct The turbulent flow in a long square duct, provides a useful validation check for the accuracy of the numerical procedures in our code, since it involves significant secondary flows, has simple boundary conditions and related experimental data have been reported. In a turbulent flow within a straight duct of rectangular cross section, a transverse mean flow exists even when the flow is fully developed. This transverse flow, commonly known as the secondary flow, brings the fluid into a lateral spiral motion superimposed upon the axial mean flow. The secondary flow of this type is termed by Prandtl as the secondary flow of the second kind. In this secondary flow the core fluid flows toward the corner along the corner bisector forming two recirculating cells when viewed in the cross-flow plane. This secondary velocity is only of the order of 2% of the streamwise bulk velocity. These secondary flows are generated due to the anisotropy of the cross-planar (transverse) normal stresses. Therefore, to predict these secondary flows more accurately, it is necessary to use a turbulence model that can account for anisotropy of the turbulence. This flow therefore is a good test case for any non-linear turbulence model. Experimental investigations on fully developed flows in rectangular ducts were carried out by Hoagland [42], Gessner [43,44], Launder and Ying [45] and Brundntt and Baines [46] . Because these particular flows provide a validation test case for turbulence models, numerical investigations have also been carried out by (among others) Launder and Ying [47], Gessner et al. [44], Nkayama et al. [48], Demuren and Rodi [49] and Myong and Kobayashi [50]. In the present study, flow simulations were carried out for a square duct with side h=0.25 m and length L=92h. Figure 2-1 shows the geometry of the duct and vector plot of the secondary flows at the end of the duct. 23 Figure 2-1: Vector plot of the secondary flows in the square duct. To insure grid independence, computations were carried out on three different grid sizes with grid 1, 2, and 3 having 8x8, 12x12 and 16x16 cells at each cross section respectively and 46 cells along the length of the duct (table 2.2). Grid 2 has approximately twice and grid 4 has four times the number of grids as the grid 1, for the same calculation domain. By taking advantage of flow symmetry along the centerlines, actual computations were carried out only for a one quadrant of the duct. 2 4 Table 2-2: Grid sizes for testing the grid independence of the square duct flow simulation Grid Cells in the X Cells in the Y Cells in the Z Total number direction Direction direction of cells 1 46 8 8 2944 2 46 12 12 6624 3 46 16 16 11776 Figure 2.2 shows the magnitude of the secondary flow along the Y axis at the X= 80 x h from the entrance for three different grids. The results in the figure 2.2 are obtained using the using the quadratic algebraic stress model described in section 2.6.1 . Figure 2.3 is similar to the figure 2.2 except that the results in the figure 2.3 are obtained using the cubic algebraic stress model described in section 2.6.2 . Figure 2.2 and 2.3 show that there is no significant difference between the numerical results obtained using the three different grid sizes and the numerical results of these flow simulations are essentially independent of the grid sizes. t i i i i i i i i i i i i i i i i ' ' ' ' 0 0.25 0.5 0.75 1 Y* Figure 2-2Secondary flow calculation using the quadratic algebraic turbulence model for three different grid sizes at X=80x h along the Y-axis for Red=250000 (D„ is the hydraulic diameter and Red is the Reynolds number based on Dh). Y* =Y/a where a is the half width of the duct. 2 5 L I I I I I I I I I I I I I I I I I I I I 0 0.25 0.5 0.75 1 Y Figure 2-3: Secondary flow calculation using the cubic algebraic turbulence model for three different grid sizes at X=80x h along the Y-axis for Red=25Q000 (D h is the hydraulic diameter and Rea is the Reynolds number based on Dh). Y* =Y/a where a is the half width of the duct. The direction and magnitude of the secondary flows along the Y axis at X=80xDh and Red=250000 (Dh is the hydraulic diameter and Red is the Reynolds number based on Dh) are compared to experimental and other numerical results available in the literature in Figure 2.4. For simplicity only the numerical results from the grid 3 are only shown in Figure 2.4. There is a good agreement between the present calculation and the numerical predictions by Myong and Kobayashi[50]. However, the experimental data reported by Gessner et al. [44] are slightly higher than the results of the present study. The largest difference between the experimental and numerical data is about 10% of the magnitude of the secondary flow itself, which is less than 2% of the bulk velocity. The secondary flow magnitude resulting from the numerical simulation using the standard k-e model is 2 6 less than 10 ^ (when it normalized with bulk velocity). This is less than the convergence criteria of the numerical simulation which required the residual error to be less than 10 " 6 . Figure 2.5 shows the magnitude of the secondary flows along the Y axis at different distance from the entrance (for simplicity only the numerical results from the quadratic nonlinear turbulence model is shown). The Figure 2.5 shows that secondary flows reach their maximum values only in the fully developed region x/d > 84. S Present Study (Launder Model) A Present Study (Spezial Model) <3 Present Study (Standard k-e model) Q Gessner and Emery (Experimental) •0 Myong And Kobayashi (Numerical) 0 0.25 0.5, 0.75 1 Y Figure 2-4 : Experimental and numerical calculation of secondary flows along the Y axis at X=80xDh and Red=250000 (Dh is the hydraulic diameter and Red is the Reynolds number based on Dh). Y* =Y/a 2 7 Figure 2-5: Secondary flow calculations along the Y axis at different distances from the entrance (quadratic nonlinear turbulence model). Y*=Y/a 2 8 3 -Experimental Arrangement and Measurement Techniques In this chapter the experimental facilities and methodology used for this study of the headbox flow are presented. The model and the water facilities where the measurements were made are described and the geometries are shown. The primary experimental measurement tools, namely the laser Doppler velocimeter (LDV) is described, along with the parameters of operation for the measurements. 3.1 Experimental Apparatus and Equipment The experiments were performed in the Pulp and Paper Center at the University of British Columbia. The experimental set-up used a closed flow system diagrammatically shown in Figure 3-1. In the flow loop, the water is pumped from the reservoir tank, which can contain a total volume of 3 m 3 of fluid, to the headbox through pipes and control valves. The flow through each tube is metered and adjusted individually. The volumetric flow rate was measured through each tube by measuring the pressure difference across Gerand engineering D7 orifice plates, placed well upstream of each diffuser tube. Experiments were conducted in a transparent plexiglass headbox, within a flow loop as shown in Figure 3-2, to allow for L D V measurements and visual inspection of the flow. This headbox was designed and built specifically for the present research and is a scaled model of a typical headbox with the size reduced by a factor of 5 in the streamwise direction and in height and the width reduced approximately by a factor of 10. The model is constructed of three different modules: tube bank, calming section and converging section connected together. This design offers more flexibility in conducting experiments for various headbox geometries (i.e. different contraction ratio and length of the calming section). In addition this design allowed for access to the inside of the headbox model for cleaning and other maintenance jobs. The rectifier tubes (or diffusers) are round at the inlet and rectangular at the outlet with slowly increasing cross sectional areas and are typically used in practice both to provide the turbulence energy and to generate a fairly uniform velocity profile at the converging section inlet. To improve the flow condition at 29 the inlet to the diffuser tubes a long straight pipe with a constant area cross section was used upstream of each diffuser. Figure 3-3 is a drawing of a typical diffuser. To avoid separation inside the diffuser tubes, the divergence angle was chosen to be around 6 degrees, below the recommended limit of 8 degrees [51]. Divergence angle is not the only parameter affecting the flow pattern inside the diffuser; the length of the diffuser also has a great effect on the type of the flow pattern inside the diffuser [51]. Altogether there are 40 rectifier tubes in the present model, two rows in the headbox height direction and 20 in the span direction. After travelling through the asymmetric converging section, the flow is finally discharged at the nozzle or "slice". The headbox section starts with a rectangular channel, which remains constant in cross sectional area until the channel length reaches 0.075 m. Downstream of this point, the channel converges to the nozzle with a total contraction ratio of 10. Details of the headbox geometry are given in Table 3-1 which gives the dimensions of the headbox converging section. One of the requirements for the dynamic similarity between the model and the prototype is the equality of the Reynolds number, which is defined as: Where U , p, p are the fluid velocity, density and viscosity respectively and h is a geometrical characteristic length. Since the model size is reduced by factor of five the velocity in the model would have to increased by factor of five to achieve the same Reynolds number in model and prototype. However the fluid velocity in the model is limited by the pump capacity. The model Reynolds number based on the average bulk velocity at the inlet (0.2 m/s), inlet height (0.075 m) and the physical properties of water at 20 °C is approximately 15000 which is much smaller than the Reynolds number of the prototype (> 500000). However, since the flow in the model and the prototype are both in the fully turbulent region the turbulence characteristics are largely self-similar and independent of the 30 Reynolds number. Therefore the flow in the model is similar to the flow in the prototype even though the Reynolds numbers are very different. To avoid entrapment of air bubbles in the model water was circulated through the headbox for at least 2 hours before taking any measurements. Table 3-1. The Geometry of the Headbox Converging Section. parameters values width 0.75 m (constant) inlet height 0.075 m slice height 0.0075 m contraction ratio 10 length 0.3 m z Figure 3-1: The Flow loop in the experiment. 31 Figure 3-2: Experimental model drawing. All dimensions are in meters. 0.001 0 . T H I C K Flow Direction 6 - ? ° 3 Z 0 . 2 3 0. J • 0 . 0 5 T S 1 — 0.0375 *— 0.0191 Flow Direction Figure 3-3 : Drawing of a typical diffuser. All dimensions are in meters. Wall thickness ~ 1 mm. 3.2 Laser Doppler Velocimetry A l l of the mean velocities and turbulence quantities for these experiments were measured using laser Doppler velocimetry (LDV). In this section the principles of operation and the particular system used for the measurements are described. 3.2.1 Background Laser Doppler velocimetry is a well-proven technique that measures fluid velocity accurately and non-invasively. Laser light illuminates the flow, and light scattered from particles in the flow is collected and processed. In practice, a single laser beam is split into two equal-intensity beams which are focused at a common point in the flow field. A n interference pattern is formed at the point where the beams intersect, defining the measuring volume. Particles moving through the measuring volume scatter light of varying intensity, some of which is collected by a photodetector. The resulting frequency measured by the 33 photodetector is related directly to particle velocity in one direction, in the plane defined by the two beams and in the direction of the bisector of the two beams. It is apparent that in this configuration it is not possible to determine the direction of travel of the particle. In order to remove this ambiguity it is possible to frequency-shift one of the beams in the pair which results in a fringe pattern which moves within the measurement volume at a speed proportional to the shift frequency. If the shift frequency is high enough, a particle passing through the sample volume will always move in the same direction relative to the fringes. The light that is scattered by a particle passing through the sample volume must be gathered by a receiver to be processed. If the receiver is incorporated into the probe that transmits the beam pair the method is referred to as 'back-scatter'. One pair of laser beams is required to measure one component of velocity. If more than one component of velocity is required, more beam pairs are needed and all must intersect at a common point. 3.2.2 Apparatus In these experiments a Coherent Innova 12.0W argon-ion laser along with a TSI Colorburst 9201 multicolour beam separator was used to generate the two, different-coloured beam pairs: green and blue. The Colorburst 9201 contains a Bragg cell to allow for frequency shifting of one of the beams in each of the beam pairs. The beams were directed through fiber optic cables to a probe, TSI model 9831, which was positioned to focus the beams to a single measurement volume. The configuration is shown in Figure 3-4 and the properties of the beam pairs are indicated in Table 3-2. The fiber optic cables allowed for flexibility in positioning the probe. The probe also contained the receiving optics, which coupled the received Doppler signal back through the fiber optic cables to a TSI Colorlink 9230 photomultiplier. The photomultiplier filtered the appropriate-coloured beams and passed the signal on to two TSI IFA550 signal processors, one for each component of velocity. 34 Table 3-2 : LDV beam properties Probe #1 Green Beam Pair Probe #1 Blue Beam Pair X [nm] 514.5 488.0 de [mm] 2.82 2.82 df [urn] 3.73 3.54 dm [Jim] 90.5 85.8 / m [mm] 1.31 1.24 K [ ° ] 3.95 3.95 where X- wavelength de - diameter of beam at probe d/ ^fringe spacing dm =diameter of the measurement volume lm - length of measurement volume K = half angle of the beam pair Beam —1 Separator —1 Laser Headbox Green/Blue Probe V Optics Beach Figure 3-4: LDV system schematics 3.2.3 Data Acquisition The acquisition and storage of the raw data was controlled by software from TSI - FIND v. 4.03. The program allows the user to select the frequency shifts, signal filters (high and low cut-off), coincidence windows, time between measurements and provides a histogram of the most recent data points. A signal was accepted as valid i f the IFA550 processors counted eight consecutive fringe crossings from a particle passing through the measurement volume, subject to some constraints on the timing of the fringe crossings. The time for the eight fringe crossings, denoted by tg; as well as the time between consecutive valid data points, tbd, were stored on disk. The velocity of the particle could then be determined from: K=<//(8/t8 x looo-yy 3.1 where V = velocity, in m/s fs = frequency shift, in MHz dj- fringe spacing, in urn tg - time for eight fringe crossings, in ns and the fringe spacing is known for each beam pair (see Table 3-2). 3.2.4 Alignment Two different types of alignment were required for these experiments: alignment of the beam pairs with respect to each other, and alignment of the beams with respect to the coordinate axis of the apparatus. Both of these issues are discussed in this section. As mentioned previously, the measurement of different velocity components requires more than one pair of beam to intersect at a common point. With the system used in these experiments the green and blue beams were transmitted through a common lens in one probe and were designed to intersect at one point with the plane of the beam pairs 36 perpendicular to each other. In this study only two components of the velocity were measured at one time, which eliminated the need for alignment of the third beam. The second alignment concern is that of the probes relative to the coordinate system of the headbox. The importance of alignment relative to the reference coordinate system may be seen by considering the measurement of a small velocity component in the presence of a large one. If the beams are not aligned accurately, part of the large velocity will be picked up by the beams which should be measuring the small component. The effect may be corrected i f the other two velocity components are known (or may be estimated) or i f the precise angle of relative alignment is known. For the complex, three-dimensional flow present in the current investigation, a priori knowledge of the flow field is not available. In this study beams are always kept at the 90° angle to the flow to minimize the error in measurements of the smaller components of velocity. Error due to misalignment angles for the present study is discussed in chapter 3.6. A detailed description of the alignment procedure, along with the measured angles of alignment, is given in Appendix A . Finally, the position of the measurement volume with respect to the channel coordinate system was determined. This alignment was performed by moving the measurement volume to several reference locations within the test section. The measurement volume was positioned at the reference location and the position of the probe traverse mechanism was recorded. 3.3 LDV operation parameters In these experiments to obtain measurements of the mean and fluctuating components of the flow field the L D V system was operated in the back - scattering mode for collecting the reflected light. Verification of the received data can be done in two ways: coincidence mode and independent mode. If the measurements of the turbulence shear stresses are necessary then the L D V system must be operated in the coincidence mode. In this mode a measurement is accepted only i f two or three IFA550 processors received a valid signal within a user-selected time span, or coincidence window. However, sometimes in this method the data rate is reduced to unacceptable levels. On the other hand, in the 37 independent mode data received from each set of beams is processed independently, increasing the data rate. In this study the L D V system was used only in the independent mode. 3.4 Seeding particles The proper selection of particles for seeding can have a strong influence on the quality of the data in a L D V system. For the forward-scatter mode in water, naturally occurring particles in the water will usually give a nearly continuous signal and particles follow the flow well because of the small size of typical particles(= 10 u\m) and the high density of the water. However in the back-scatter mode similar to the one used in the present study seeding particles were needed to improve the signal to noise ratio and the data rate. The parameters that determine whether the particle is following the flow are the size, density and shape of the particle. Likewise, the only other parameter of concern in the velocity measurement is the strength of the scattered light or, more precisely, the signal to noise ratio of the resulting signal from the photodetector. Acceptable particle diameters for the amplitude rjs | U p / Uf | =0.99 (where the Up and Uf are the instantaneous velocities of the particle and the fluid respectively) is dependent on the density ratio of the particle to the density of the fluid. When the ratio of the particle density to fluid density is equal to unity, particles will follow the flow exactly regardless of the size (aside from the effect on the flow of the particles themselves). In the present study to reduce the noise due to the reflection from the walls, seeding particles with a high refractive index were needed. For this reason metallic coated hollow glass spheres were used to seed the water for the L D V measurements. These coated hollow glass spheres were selected for seeding due to their low density and high reflectivity. The properties of these particles are given in Table 3-3. A n acceptable particle diameter in water for the density ratio of the particles used in this experiment is between 5 |im and 16 \im Durst [52]. 38 Table 3-3: Properties of the metallic coated hollow glass spheres. Mean diameter(um) Size range(um) shape Specific gravity 10 10% <7 90%<21 Spherical 1.65 In addition to the properties of the seeding particles, other factors must be taken into account for the L D V measurements performed in these experiments. The first involves possible errors from non-uniform seeding of the flow. One type of non-uniform seeding occurs when the seeding density is dependent on the local flow velocity. There may be regions of densely-seeded, low velocity flow, which would result in more low velocity readings being made by the processors, thus biasing the ensemble averages used to calculate the velocities. This form of non-uniform seeding was eliminated by thoroughly mixing the seed particles with the flow before entering the model. The second factor which must be taken into account is the effect of multiple realizations. The IFA550 processors continuously check the incoming signal to determine whether or not a valid Doppler signal is present. After each measurement of tg is accepted, the processor resets to accept another signal, but i f a seed particle remains inside the measurement volume during this time the particle velocity may be measured again. Slower-moving particles reside in the measurement volume for a longer time and will produce a greater number of realizations (in the limit of a zero velocity particle where frequency shifting is used the number of realizations is infinite) which will bias the flow statistics. This can be shown by considering an example U=1.0m/s .£=100 kHz df= 3.73 urn <im=90.5 um where the U is the component of the velocity in the direction measured by the green beam pair. Using the equation 3.1 we can calculate t8=21734 ns. If the particle travels the full diameter of the measurement volume it will have a residence time of 90500 ns (dm/U) which could possibly result in 4 readings(=90500/21734) from the same particle. 39 The selection of a measurement window which was close to the expected residence time(calculated as above) forced the processors to reset after the measurement window time, rather than immediately after tg was measured. This helped to ensure that a given particle was no longer in the measurement volume. The bias in the mean flow data due to fluctuations in the velocity, which is often discussed in the L D V literature, is discussed in section 3.5.2 of this thesis. 3.5 LDV Data Analysis As discussed in section 3.2.3, the data from each measurement position in the flow field was stored on disk in the form of t» values for each data point. At each measurement point data were taken for 30 s or until at least 800 data point were obtained. These values were converted into velocities using the FIND 4.3 software program [53]. The analysis consists of the determination of the average velocity, rms velocity, skewness coefficient and flatness coefficient. 3.5.1 Removal of Bad Data Points The FIND software allows the user to remove velocities which are outside three standard deviations from the mean. This allowed the analysis to exclude data points which were considered strange or unlikely and was used for all the results presented here. The removal of data outside three standard deviations has been used previously [54] to remove the apparent turbulence produced by the ambiguity noise in the signal [55]. Tests for the experiments performed here indicated that there was no noticeable difference between the mean velocity or the rms values of the fluctuating components of the velocity results obtained using this approach and that with no data point removal. 3.5.2 Calculation of Expected Values - Velocity Bias Correction When mean flow statistics are calculated from L D V data using arithmetic averages, a bias in the statistics will exist [56]. The velocity bias described here differs from the bias 40 described in section 3.4 in that it is assumed the flow is uniformly seeded. The bias described here occurs because, in a uniformly seeded flow, more particles will pass through the measurement volume in a given time when the flow velocity is high than when the velocity is low. If a simple arithmetic average is taken, the statistics will be biased towards the high velocities. A number of different correction schemes, typically weighted averages, are available to reduce this bias [56]. FIND software uses transit time weighting to correct for velocity bias. In the transit time weighting the total time burst x is also recorded along with the velocity of the particle. In the IFA 550 the total time burst is equal to the time required for the particle to cross the 8 fringes tg, in the equation 3.1. Then T is used for statistical processing of the data as follow; Velocity mean (f)=-Standard deviation {a 1 = Turbulence intensity (%) T% = It-o\,xl00 a detailed description of the method can be found in the TSI FIND 4.3 software manual [53]. Even with a good velocity bias correction scheme such as transit-time which was used in the present study errors may occur due to the turbulent nature of the flow [57]. In their simulations of a turbulent flow field the percent error in the mean velocity and variance is defined as: (lj)—U (u'u')-uu' = ^ - L and 62-U w V where <) indicates the estimate of the true value. When ensemble averaging was used, it was found that pi increased with the square of the turbulence intensity (T 2 = uu 2/U 2) for turbulence intensities below 30%. This value represents the maximum likely error in the 41 mean velocity. When transit-time weighting was applied the error (3 i was greatly reduced (pV 5%) but p 2 could still be relatively large ( p 2 « 20%). While differences exist between the parameters used in the simulations of Fuchs et al [57] and those present in the current experiments, the results do provide a sense of the magnitude of the bias caused by turbulence and the level to which the bias may be corrected. Transit time weighting, which was used in the present study, is a reasonable method to use for dealing with velocity bias but errors will still remain, with the magnitude varying with the local flow conditions. However, additional corrections beyond the initial weighting scheme are difficult to justify and the errors induced by the turbulent flow are probably the primary cause of uncertainty in the experiments. 42 4 - Experimental measurements In this section L D V measurements of mean and fluctuating components of the velocity in the headbox model described in the chapter 3.1 of this thesis will be presented and discussed. The techniques described in the chapter 3 of this thesis are applied for collecting and processing the data at each measurement point. 4.1 Error analysis In this section we present a discussion of the types of error that may be present in the experimental data. The real errors in experimental data are those factors that are always vague to some extent and carry some amount of uncertainty. Therefore we need to determine just how uncertain a particular observation may be and to show this uncertainty in the experimental data in a consistent way. There are two types of error that may cause uncertainty in the experimental measurements, fixed or systematic errors and random errors. Systematic errors will cause repeated readings to be in error by roughly the same amount but for some unknown reason. On the other hand random errors may be caused by personal bias, random electronic fluctuations in the apparatus or instruments, etc. These random error usually follow a certain statistical distribution but not always. In general in most instances it is very difficult to distinguish between fixed errors and random errors. Some general rules can be used to get a first estimate of the uncertainty in the final results due to the uncertainties in the primary measurements. For example one can assume the error in the result is equal to the maximum error in any parameter used to calculate the result. Kline [58] has presented a more precise method of estimating uncertainty in experimental results. The method is based on a careful specification of the uncertainties in the various primary experimental measurements. Suppose a set of measurements is made and the uncertainty in each measurement can be expressed. These measurements are then used to calculate some desired result of the experiments. We wish to estimate the uncertainty in the calculated result on the basis of the uncertainties in the primary 43 measurements. The result R is a given function of the independent variables x i , X2 , x n . Thus, R=R(xi ,x 2 , ... x n) Let w r be the uncertainty in the result and w J ; w2, , w n be the uncertainty .in the independent variables. Then the uncertainty in the result is given as[58] Wr=\ r dR \2 f w l J dR w2 J Kdxn \ w n J 4.1 In the following sections, the error in each primary measurements affecting the L D V measurements will be discussed. 4.1.1 Angular alignment of the LDV system As was discussed in section 3.2.4, the probe was aligned at a 90 degree angle relative to the headbox coordinate system. The angle of alignment was determined using the techniques described in Appendix A . The variance of each measured angle is used to estimate the error in the velocities. More details on the approach used here for error calculation due to angular alignment may be found in Findlay [59]. As an overall estimate of the error, the mean velocities are typically ± 0 . 0 1 U (mean streamwise component of the velocity) and the uncertainties in the turbulence normal stresses are +0.001 U (mean streamwise component of the velocity). 4.1.2 Measurement volume positions The position of the measurement volume in the model coordinate system was determined after the probe was aligned (see section 3.2.4). The standard deviation of the repeated measurements gives the positioning accuracy of the measurement volume relative to the model coordinate system. The largest value obtained for the positioning error was 0.7 mm which is 1.8% of the tube width at the exit. 44 4.1.3 Bulk velocity measurements In this experiment the bulk velocity through each diffuser tube is measured and controlled separately by measuring the volumetric flow rate as described in section 3.1. Volumetric flow rate is measured using a differential pressure meter to measure pressure across an orifice plate. The flow meters were calibrated by collecting and measuring the volume of the water flow through flow meters for 3 minutes using a stop-watch and a measuring cylinder. Then this volume flow rate was compared to the value obtained by converting the pressure drop across the flow meter to volume flow rate using the chart provided by the company. Although the manufacturer of the system claims the accuracy of ± 1%, the random testing of the flow meter showed the error in volumetric flow rate measurements of approximately + 4%. The other factor affecting the bulk velocity measurement is the exit area of the diffuser. Measurements of the exit areas of several diffusers showed the mean value of area as 1406.25 mm 2 and standard deviation of 36.5 mm 2. Using equation 4.1 we can calculate the uncertainty in the measurement of the bulk velocity as follow: U= Volume flow rate(Vf) / Area 3U/3Vf= 1/Area 3U/3Area= Vf/Area 2 Area=0.00140625 ± 0.0000365 m 2 V f 0.0003375 ± 4% m3/s w A= 0.0000365 m 2 wv=0.0000135 m3/s W u={((8U/9V f )x w v ) 2 + ((3U/aArea)xwA)2 } 1 / 2 W u = {(l/0.00140625)x0.0000135)2+((0.0003375/ 0.001406252 )x 0.0000365)2 } 1 / 2 W u = 0.01144 m/s This is + 5 % of the theoretical bulk velocity of the 0.24 m/s. Therefore the variation of the bulk velocity at the exit plane of the different diffusers is ± 5.0%. Figure 4.1 shows the streamwise component of the mean velocity measured horizontally (in the cross machine direction) at the centerline of the four diffusers on the bottom row of the 45 diffusers. L * is the distance in the cross machine direction measured from the centerline on the center diffuser divided by the width of a diffuser. Figure 4-1 shows the exit velocity through each diffuser. Bulk velocity through each diffuser in Figure 4-1 are tabulated in Table 4-1. Comparison of the values in the table shows the difference between the bulk velocity through each diffuser falls with in the expected experimental error. Values in the table are obtained by integrating the velocity profiles for each diffuser in the Figure 4-1 using the trapezoid integration method. Figure 4-1: Mean velocity profile measurement across different diffusers. Data are measured at the centerline of the diffusers on the bottom row. L * is the distance from the centerline of the center diffuser at the bottom row divided by the width of a diffuser. Table 4-1: Bulk velocity through different diffusers at the bottom row. Tube number Bulk velocity (m/s) 1 0.2573000 2 0.2500834 3 0.2568500 4 0.2619335 46 4.2 Velocity measurement results The length and the divergence angle of the diffusers for this experimental apparatus were chosen to provide a steady viscous flow similar to that in the Figure (4.2 a). The measured flow patterns were not like those of Figure (4.2 a) or Figure (4.2 b) suggesting that there may have been some stall and transitory flow with in the tubes. Manufacturing defects may be the cause of this unexpected flow pattern leaving the tubes. A measurement of the streamwise component of the flow 3-cm after the diffuser exit is shown in Figure 4-3. In Figure 4-3, L * is the distance from the diffuser side wall divided by the width of the diffuser. Figure 4-4 shows measurement of the velocity profile done vertically at the centerline of a diffuser column. Although the volume flow rates through each diffuser is the same, for the reasons not known to the author at the present, the velocity profile is different for the top and the bottom diffuser. This velocity pattern was observed in all the diffusers. To ensure that this behavior is not due to the difference in the flow path leading to the diffusers the pipes for the top and the bottom diffusers where changed but the difference in the velocity profiles remained the same. This unexpected pattern needs further investigation but is unlikely to significantly affect the overall flow in the contracting section. Figure 4-6 shows the mean components of the velocity measured in the region after the diffusers at different locations along the center line of the symmetry plane (line A B C in Figure 4-5). Figure 4-7 shows measured root mean square values of the fluctuating components of the velocity in the X , Y and Z direction u ' , v ' and w' normalized with the local mean component of the velocity in the X direction (U). In the Figure 4-8 the rms data are normalized by the value of u ' at the inlet (u'i n). Rms values of the fluctuating component of the velocity in the X direction, u 1 diminish towards the exit, while w ' increases. Fluctuating components of the velocity in the Y direction v' remains approximately constant. Figure 4-8 shows that turbulence field becomes roughly isotropic a short distance inside the converging section. The turbulence field becomes very non-isotropic as the flow nears the exit. 47 point Figure 4-2 : Flow pattern inside a diffuser ( from White [51]). 0.3 0.25 0.2, w E =0.15 0.1 0.05 5 I _l I I L 0.25 0.5 L* 0.75 _ i i i i _ Figure 4-3: Flow measurements across the width of the diffuser. Measurements are done 30 mm downstream of the exit plane of the diffuser. L * is the distance from the diffuser side wall nondimensionalized with the width of the diffuser. U is the streamwise component of the velocity. 0.8 h 0.6 h 0.4 h 0.2 0 hEH I'I—11 n=ni r - i _ i i i i L 1 1 1 u 0 0.1 0.2 0.3 U m/s Figure 4-4: Streamwise velocity profile measurements at the centerline of a diffuser column. Measurements done 30 mm downstream of diffuser exit plane. Z* is the height from the bottom plane of the headbox divided by the total height of the headbox at that position. 49 Turbulence models which are based on the assumption of isotropy, are not likely to give accurate estimates of the fluctuating components of the velocity or the overall kinetic energy, specially in the region close to the exit. Velocity measurements were continued for a small distance into the free surface jet. Figure 4-9 shows again the variation of the mean streamwise component of the velocity in the mid plane along the line A B C but in the Figure 4-9 the measurements were continued for a short distance into the free surface jet. Experimental measurements in the free surface jet are done only for the short distance after the exit plane for two components of the velocity namely U (streamwise direction) and V (cross machine direction). Due to the short thickness of the jet and presence of surface waves accurate measurements of the W component of the velocity were not possible. In the Figure 4-9 U * is the mean streamwise component of the velocity divided by the bulk velocity at the inlet to the headbox section and X * is the distance along the X-axis divided by the total length of the converging section. Figure 4-9 shows that the U * values do not decrease immediately after the exit plane indicating the existence of a short contraction after the exit. Figure 4-10 shows the root mean square values (rms) of the velocity fluctuations normalized by the rms value of the velocity fluctuation in the X direction (u 1 m ) at the inlet at the same positions as Figure 4-9. Figure 4-5: Measurement location along the center line of the symmetry plane. 50 10 8 * 6 tL 4 2 0 -2 0 -0 Measurements W* - B Measurements U* H> Measurements V* J I I L 0.25 J l _ L j L 0.5 0.75 Figure 4-6: Mean velocity components measured experimentally in the region downstream the diffusers, along the line ABC in figure 4.5. Velocities are nondimensionalized with the bulk velocity at the inlet to the converging section. X * =X/Xi where X i is the total length after the diffuser tubes. 51 0.3 -0.25 -| 0.2 h 0.1 -0.05 -6 1 1 1 1 • • 6oc o 6o • i 0.2 1 1 1 1 0.4 _ i i i >_ 0.6 Q<X) D 6 J i i i 0.8 Figure 4-7: Rms. values of the fluctuating component of velocities u', v', normalized by U (local mean velocity in the X direction) along the line ABC figure 4.5. 1.6 Figure 4-8: Experimental measurements of rms. values of the fluctuating component of velocities u' , v' and w' normalized by u 1 inlet along the line ABC in figure 4.5. 5 3 Figure 4-9: Mean component of the velocity in the streamwise direction nondimensionalized with the inlet bulk velocity. X* is the distance along the X axis normalized with the total length of the converging section 5 4 1.2 1 0.8 0.6 0.4 0.2 0 O • Measurements u'rms o Measurements v'rms • o o rt o e o O G o 0.2 • • Exit Plane 0.4 i 1 1 — 0.6 • ' ' 1 1 0.8 1 X* Figure 4-10: Rms. values of the fluctuating component of velocities u", v', normalized by u' inlet. In addition to the measurements along the center line, velocity profiles were measured along the vertical lines at different distances along the X axis in the symmetry plane. Figure 4-11 shows the schematic of these measurements locations. Figure 4-12 shows the mean component of velocity in the X direction along the vertical line at different positions along the X axis shown in the Figure 4-11. The Z* in the Figure 4-12 is the distance from the bottom plane normalized with the total height at the same location. Figure 4-12 shows that the velocity profile variation in the Z direction (paper thickness direction) reduces significantly towards the exit and the foot prints of the individual diffusers vanish within a short distance into the converging section. Closer examination of the velocity profiles along the vertical lines are shown in Figure 4-13 which shows that although the footprints of the tubes disappears within a short distance, there is a slight variation of velocity in the Z direction even close to the exit. Fluctuating components of the velocities in the X direction (u 1), Y direction (v ' ) and in the Z direction(w') at same measurement points are shown in Figure 4-14 , Figure 4-15 and Figure 4-16 respectively. These figures show that after a short distance into the converging section, variation of the fluctuating components of the velocities along the vertical lines reduces significantly except for small regions close to the walls. 55 B Figure 4-11: Schematic drawing of the velocity profile measurements location. Measurements are done along the vertical line at the center of a diffuser set (line AB) in the X direction. 56 1 r 0.9 -0.8 -0.7 -0.6 -N 0.5 r 0.4 -0.3 -0.2 -0.1 -0 0 • A o o X= 3 cm X= 9.4 cm X= 17.5 cm X= 27.5 cm 0.2 0.4 U(m/s) 0.6 O o o o 0.8 Figure 4-12: Mean velocity component in the X direction measured along the line AB in figure 4.9 at different X position. Z* is the distance from the bottom plane normalized with total height at the same location. 57 ir-0.8 I-0.7 I-0.6 r 0.4 I-0.3^ 0.2 I-0.1 r X= 3 cm • • • • • 0.05 0.1 0.15 0.2 U(m/s) 1r 0.9 0.8 0.7 0.61 N0.5 0.4 0.3 0.2 0.1 0 • • X= 9.4 cm • i . . i i , , , , i , , , , i 0.19 0.2 0.21 0.22 0.23 0.24 U(m/s) 1 0.9 0.8 0.7 0.6 N0.5 0.4 0.3 0.2 0.1 0 • X= 17.5 cm • • 0.305 0.31 0.315 0.32 U(m/s) 1 0.9 0.8 0.7 0.6 N0.5 0.4 0.3 0.2 0.1 0 z X= 27.5 cm 0.76 0.765 0.77 U(m/s) _i_ 0.775 0.78 Figure 4-13: Variation of the velocity profile in the Z direction at different location along the X axis. Z* is the distance from the bottom plane normalized with total height at the same location. Measurements are along the line AB in figure 4.9. 5 8 1 0.9 0.8 0.7 0.6 0.5 0.4 IT 0.3 -0.2 -0.1 -0 o o o o o o o o o o o o o o o o o A A A A A • A A A A A A A A A • • • • • • • • • • • • • • • • • • • 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 iT rms(m/s) X= 3 cm X= 9.4 cm X= 17.5 cm X= 27.5 cm Figure 4-14: Fluctuating component of velocity in the X direction (MD) measured along the line AB in figure 4.9 at different X position. Z* is the distance from the bottom plane normalized with total height at the same location. 59 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 .o o o ; _ • X= 3 cm ; A X= 9.4 cm Z 0 X= 17.5 cm o X= 27.5 cm o o o o o • o A A A A A A A A A A A A A A A j i i i _ _l I I L . A A A A i i • • • • • • • • • • • • • • • • • 0 0.01 0.02 0.03 v" rms (m/s) 0.04 0.05 Figure 4-15: Fluctuating component of velocity in the Y direction (CD) measured along the line AB in figure 4.9 at different X position. Z* is the distance from the bottom plane normalized with total height at the same location. 60 0.9 0.8 0.7 0.6 N 0.5 0.4 0.3 0.2 0.1 0 • A o o X= 3 cm X= 9.4 cm X= 17.5 cm X= 27.5 cm A A A 8 A A A A A A A A A C • • • • O • 0 o • o D o • On 0.01 0.02 n 0.03 ^ 0.04 w* rms (m/s) O o o 0.05 0.06 - i i i i_ Figure 4-16: Fluctuating component of velocity in the Z direction measured along the line ABC in figure 4.5 at different X positions. Z* is the distance from the bottom plane normalized with total height at the same location. 6 1 5 - Numerical Simulations 5.1 Turbulent flow through the converging section of a headbox This section of the thesis focuses on numerical simulations of the mean and turbulent flow through the converging section of a headbox and their comparison with present measurements of these quantities. The effects of the different boundary conditions used for numerical simulations are discussed in chapter 5.1.1. Both the standard k-e and non-linear k-e turbulence models described in previous chapters were used to simulate the turbulence generated secondary flows in the converging section of a headbox as described in chapter 5.1.3. A first simulation was made to compare with the previous numerical results published by Aidun and Kovacs [20] The geometry and boundary conditions for the numerical simulations are the same as those used by Aidun and Kovacs [20]. Figure 5-1 and Table 2-1 describe the computational geometry. Table 5-1: Computational geometry basic dimensions. Length in X direction Lx=1.125m Width in Y direction LY=3.75 m Height at inlet in Z direction Lz=0.375 m Contraction ratio 10 62 ' 0.0375 Figure 5-1: Computational geometry. All dimensions are in meters. 5.1.1 Effects of the different boundary conditions The choice of boundaries requires special attention. The exact boundary conditions are not known at the exit. There are a few possible inaccurate boundary conditions that can be chosen and evaluated, one being a constant uniform velocity. This boundary condition is clearly not accurate and is not recommended. The second boundary condition is the modified gradient boundary condition, which allows for the slice contraction. In this modified gradient boundary condition, gradients of the streamwise component of the velocity U for each computational cell at the boundary is calculated as follow dU V dA dX a dX Where the A is the projection area of the cell face in the streamwise direction. In this case dA/dX can be viewed as the contraction ratio of the cell in the X direction. For the scalar equations the gradient in the X direction is set equal to zero. dX where § is any scalar property of the fluid, such as k,e. 63 The third boundary condition is a static pressure boundary. Although the pressure condition is probably more realistic, it is more difficult numerically to converge than the gradient boundary condition. To examine the effect of these last two different boundary conditions at the exit, both pressure and gradient boundary conditions were applied to the exit of the slice with the uniform inlet boundary condition such that For the plane X=0 , 0<Y<3.75 and 0<Z>0.375 the U=constant=1.22 m/s, k=constant=0.022 m2/s2 e=constant=0.005 m2/s3 while the other conditions were kept the same. This test was repeated for the case of the non-uniform boundary condition at the inlet. In this case flow in a small region near the wall at the inlet plane was set to zero such that For the plane X=0 , 0<Y<3.7m and 0<Z>0.375 the U=constant=1.22 m/s, k= constant =0.022 m2/s2 e= constant =0.005 m2/s3 and for the plane X=0 , 3.7m<Y<3.75m and 0<Z>0.375m the wall boundary condition (described in section 2.5.1) was applied. The inlet geometry for the case of the non-uniform inlet flow is shown in Figure 5-2. This non-uniform inlet condition is selected as a simple representation of a non-uniform flow in the real headbox, which usually originates from the manifold. The geometry of the computational domain is same as the one shown in figure 5.1. Figure 5.3 shows the numerical grids in the XZ plane. Table 5.2 shows the boundary conditions applied at the other boundaries for these tests. The computations have been carried out using the CMGFD code with the standard k-e turbulence model. The results obtained for the case of the uniform inlet boundary condition are shown in the Figure 5-4, 5.5 and 5.6. These results were obtained first by applying a modified gradient boundary condition at the exit and secondly a pressure boundary conditions. 64 Table 5-2: Boundary condition for numerical calculations. Exit plane {at X=1.125 m,0<Y<3.75m and 0<Z<0.0375m } Pressure Boundary condition P e xiF 101 kPa , Or modified gradient boundary conditions Top plane {0<X< 1.125, 0<Y<3.75m and Z=0.375-Xxtanoc where a is the convergence angle=16.7 degree} Wall boundary condition Bottom plane {Z=0,0<X<1.125, 0<Y<3.75m} Wall boundary condition (described in section 2.5.1) South {at Y=0.0,(KX<1.125, 0<Z<(0.375-(Xx tan a))} where a is the convergence angle=16.7 degree Symmetry boundary condition, dV/dY =0 North {at Y=LY=3.75m ,0<X<1.125, 0<Z<(0.375-(Xx tan a)) where a is the convergence angle=16.7 degree} Wall boundary condition Figure 5.4, 5.5 show the mean components of the velocity in the X and Z direction respectively in the symmetry plane (Y=0) along the line joining the center point of the inlet to the center point of the exit. The Y component (cross machine direction (CD) of the velocity in the symmetry plane is zero and is not shown in the figures. Figure 5.6 and 5.7 shows the variation of the mean component of the velocity in the Y direction (CD direction) along the center line of the exit plane (X=1.125) for the case of the uniform inlet condition and non-uniform inlet condition respectively. One can see that in the case of uniform inlet boundary condition the two different boundary conditions at the exit of the slice give essentially the same results. The explanation is that the flow is governed almost entirely by upstream conditions and therefore the downstream conditions have very limited influence. However Figure 5-7 shows that these two boundary conditions (gradient boundary condition and pressure boundary condition) would produce different results if the inlet boundary conditions were not uniform. Therefore the use of the gradient boundary condition should be limited to the case of the uniform inlet boundary condition. 65 2r> 1.8 1.6 1.4 _ 1.2 £ 1 N 0.8 0.6 0.4 0.2 p 0, Blocked Area 2 Y(m) Figure 5-2: Geometry of the inlet for the case of the non-uniform inlet flow. 0.6 -0.5 r 0.4 -Figure 5-3 : Computational mesh at the symmetry plane. 6 6 12r Figure 5-4: Mean streamwise component of the velocity at the symmetry plane Y=0 for two different boundary condition at the exit with the uniform flow at the inlet along the line joining the center of the inlet plane to the center of the exit plane. U* is the mean streamwise component of the velocity divided by the bulk inlet velocity. X* is the distance in the X direction divided by the total length in the X direction. X* Figure 5-5: Mean component of the velocity in the Z direction at the symmetry plane Y=0. for two different boundary condition at the exit with the uniform flow at the inlet along the line joining the center of the inlet plane to the center of the exit plane. W* is the component of the velocity in the Z direction divided by the bulk inlet velocity. X* is the distance in the X direction divided by the total length in the X direction. 67 01 -0.002 h > -0.006 h -0.004 h e e Pressure Boundary Condition Gradient Boundary Condition -0.008 -0.01 0 0.25 0.5 Y* 0.75 Figure 5-6: Mean component of the velocity in the Y direction at the center line of the exit plane for two different boundary condition at the exit with the uniform flow at the inlet. V* is the mean component of the velocity in the Y direction divided by the bulk inlet velocity. Y* is the distance in the Y direction divided by the total length in the Y direction. 68 z B -0.0025 E-* >-0.005 -0.0075 F--0.01 --0.0125 b -0.015 0 - B P ressure Boundary Cond i t ion -0 Grad ient Boundary Cond i t ion 0.25 0.5 Y * 0.75 Figure 5-7: Mean component of the velocity in the Y direction at the center line of the exit plane (Line AB above) for the case of the non-uniform flow at the inlet. V* is the mean component of the velocity in the Y direction divided by the bulk inlet velocity. Y* is the distance in the Y direction from the center of the line AB divided by the total length in the Y direction. Both symmetric and non-symmetric converging sections are in use in industrial headboxes. To examine the effects of these two different geometries on the flow pattern and the generation of the secondary flow, numerical simulations of the flow in both symmetric and non-symmetric converging sections were carried out here. The geometry and the boundary conditions used for the numerical simulations are the same as those described in section 5.1. The symmetric converging section has the same dimensions and same convergence ratio as the non-symmetric one except that in the symmetric converging case both the top and bottom walls have an equal angle with 69 respect to the axis of symmetry. Simulations were done using the standard k-e model. A uniform velocity boundary condition was used at the inlet and the pressure boundary condition was used at the exit. Figure 5-8 and Figure 5-9 show the contour plot of the mean component of the velocity in the Y direction (CD direction) for the non symmetric and symmetric convergence respectively. The pattern and magnitude of the velocity is very similar for both the symmetric and non-symmetric converging geometry. However for the case of non uniform inlet condition where a small region close to the wall at the inlet is blocked (see Figure 5.2) these two geometries produce different flow patterns with the maximum value of the cross flow component of the velocity being slightly lower in the symmetric convergence case. In the case of the non-uniform inlet a cork screw motion is created close to the wall, which covers the whole height of the headbox. This cork screw motion is shown in the Figure 5.11. This cork screw motion was not present in the symmetric converging section with the non-uniform inlet. Figures 5.10 and 5.12 show the contour plot of the mean component of the velocity in the Y direction (CD direction) with the non uniform inlet condition for the non symmetric and symmetric convergence respectively. The conclusion from the numerical simulation is that non- uniform inlet conditions may produce unexpected and extensive secondary flows at the exit of a converging section. 70 Figure 5-8 : Cross flow contour plot for the non symmetric converging section with uniform inlet condition. 7 1 7 2 7 3 Figure 5-11: Cork screw motion in the non-symmetric convergence V 0.101371 0.0920412 0.0827119 0.0733825 0.0640531 0.0547237 0.0453944 0.036065 0.0267356 0.0174062 0.00807687 0.0012525 0.0105819 0.0199113 0.0292406 Figure 5-12: Cross flow contour plot for the symmetric converging section with non-uniform inlet condition. 75 5.1.2 Grid independence To insure the numerical calculations presented in this section are independent of the grid size the numerical simulations are carried out on three different grid sizes. The number of cells in each coordinate direction is given in table 5.3. Grid stretching techniques were used to concentrate the grid points in the vicinity of the side wall and toward the exit. Grid set 3 has eight times and grid set 2 has four times the number of grid points in grid set 1. Table 5-3: Grid sizes for testing the grid independence of the flow simulation in the converging section of a headbox Grid Number of cells Number of cells Number of cells total number of in X direction in Y direction in Z direction grids 1 9 30 9 2430 2 15 50 15 11250 3 18 60 18 19440 Grid sensitivity tests were carried out for all the turbulence models described in chapter 2, namely the standard k-e model, quadratic nonlinear k-e model, cubic nonlinear k-e model and Reynolds stress model. Figure 5.13 shows a comparison of the results obtained using the different grids for the calculation of the mean component of velocity in the X direction. Data are selected along the line joining the center of the inlet plane to the center of the exit plane at the spanwise symmetry plane. Figure 5.14 shows the turbulence kinetic energy normalized by its inlet values at the same locations as figure 5.13. For simplicity only the results from the Reynolds stress models are shown. Figure 5.15 shows the mean component of the velocity in the Y direction at the centerline of the exit plane. There is a significant difference between the results of the grid 1 and those from grid 2 and 3. This indicates 76 that the results using the grid 1 are not reliable and only the solution for the grids 2 and 3 are grid independent. Therefore in this section only the results from the grid 3 are presented. 10 9 8 7 6 r 5 4 t-3 \ 2 r> 11 0 0 • Grid 9 x 30 x 9 Grid 15 x 50 x 15 Grid 18x60x18 0.25 o 0.5 X* 0.75 Figure 5-13: Comparison of the mean component of velocity at the centerline for different grids. U*=U/Uin, X*=X/LX (see table (4.3 and 4.1) for definition of terms). Results obtained using RSM. 7 7 15 10 • Grid 9 x 30 x 9 Grid 15x50x15 O Grid 18x60x18 3? 0.25 _l_ 0.5 X* 0.75 • • 9 o 1 Figure 5-14: Comparison of the turbulence kinetic energy, k* = k/kin. X*=X/Lx, Here kj. is the turbulent kinetic energy at the inlet to the headbox. Results obtained using RSM. Figure 5-15: Comparison of the CD velocity at the exit plane, V*=V/Uj„ , Y*=Y/L Z, Uj„ is the bulk velocity in the streamwise direction at the inlet. Results obtained using RSM. 7 8 5.1.3 Secondary flow calculations Numerical simulation of the flow in the headbox is presented in this section. Different turbulence models were used and the results were compared. The geometry and the dimensions of the computational domains are the same as the one described in section 5.1. A uniform velocity boundary condition is used at the inlet to study the effects of the turbulence generated secondary flows isolated from any upstream effects. Although this is different from reality it is necessary to insure that secondary flows in the jet are generated in the headbox and are not due to non-uniformity of conditions at the inlet. A uniform boundary condition at the inlet is also used for both the turbulence kinetic energy and its rate of dissipation. The modified gradient boundary condition described in section 5.1.1 is applied at the exit plane for easy convergence. Boundary conditions used for these numerical simulations are summarized in Table 5-4. The mean component of the velocity in the X direction (U) is shown in Figure 5-16. Data are selected along the line joining the center of the inlet plane to the center of the exit plane at the spanwise symmetry plane. All turbulence models predict essentially the same mean velocity distribution. The mean velocity component in the Y direction (V or cross machine direction velocity) at the exit plane is shown in Figure 5-17. Results from standard k-e model and the nonlinear model are very close to each other. If these CD (y direction) velocities were due to turbulence anisotropy, and have the same origin as the secondary flows in the square duct, then one can assume they should have the same pattern and magnitude. That is, fluid should move toward the corner along the corner bisector and move toward the center along the side wall bisector. Thus, the CD velocity along the wall bisector should always be negative. Also, secondary flows in a square duct reach their maximum value of approximately 2 percent of the bulk velocity only in the fully developed region and are much smaller in magnitude in the developing region. In view of these observations from the channel flow simulation, the CD velocity magnitude and direction calculated by the nonlinear model appears to be most accurate. In this study the geometry and boundary conditions have been selected to be the same as those used by Aidun et al [20], as far as possible, but there is uncertainty regarding the exact inlet condition used by Aidun et al 79 [20]. This may be the cause of the considerable difference between the present results and those reported by Aidun et al [20]. Table 5-4: Boundary condition for numerical calculations. Inlet plane { X=0 , 0<Y<3.75 and 0<Z>0.375} U=constant=1.22 m/s , k= constant =0.022 m2/s2 e= constant =0.005 m2/s3 Exit plane {at X=1.125 m,0<Y<3.75m and 0<Z<0.0375m } Modified gradient boundary conditions Top plane {0<X<1.125, 0<Y<3.75m and Z=0.375-Xxtana where a is the convergence angle=16.7 degree} Wall boundary condition (described in section 2.5.1) Bottom plane {Z=0,0<X<1.125, 0<Y<3.75m} Wall boundary condition South {at Y=0.0,0<X<1.125, 0<Z<(0.375-(Xx tan a))} where a is the convergence angle=16.7 degree Symmetry boundary condition, dV/dY =0 North {at Y=LY=3.75m ,0<X<1.125, 0<Z<(0.375-(Xx tan a)) where a is the convergence angle=16.7 degree} Wall boundary condition 80 Figure 5-16: Mean component of velocity at the centerline. U*=U/Uj„, X*=X/LX (see table (4.3 and 4.1) for definition of terms). 8 1 L—i—i—i—i—l—i—i i i l i i i i I i i i i I i i i i I 0 2 4 _ 6 8 10 Y Figure 5-17: CD velocity at the exit plane, V*=V/U i n , Y*=Y/LZ, U i n is the bulk velocity in the streamwise direction at the inlet. 8 2 Figure 5-18: Turbulence Kinetic energy, k* = k/kj„. X*=X/Lx , Here kj„ is the turbulent kinetic energy at the inlet to the headbox. In Figure 5-18 the turbulence kinetic energy k at the symmetry plane, calculated by the different turbulence models, is shown. One has to consider that the turbulence level predicted by the numerical methods depends on the inlet e-values, which have not yet been specified accurately for the numerical simulations. Figure 5-18 is based on an inlet value of e=0.005 m 2/s 3. (This will be specified more accurately to match observed turbulence inlet values in chapter 5.2). However in general, the values from both the standard k-e and the nonlinear k-e are much larger than to the R S M values. In the case of these flows with such a high rate of longitudinal strain, the generation term in the k equation is apparently modeled in a way that produces results for k which are much too large. This problem will be addressed again in the chapter 5.3. Comparison of the experimental and numerical data shows that both the standard k-e and the non-linear k-e model predict the components of the mean velocity reasonably well. Therefore, as far as mean velocity calculations are concerned, the standard k-e model can be used for simulation of the industrial scale headboxes and manifolds. But the accurate 8 3 simulation of turbulence will require modifications of the usual engineering turbulence models. 5.2-Flow through the headbox model In this section, the numerical simulations are done on a domain similar to those of the experimental model described in the chapter 3. Figure 5-19 shows the computational domain. The dimensions of the computational domains are the same as the dimensions of the experimental apparatus shown in figure 3.2 of chapter 3. To reduce the size of the computational domain, only one column of the diffuser tubes is simulated. A uniform diffuser wall with a thickness of 1.0 mm is considered for this calculation, consistent with the experimental arrangements. Figure 5-20 shows the computational grid in the XZ cross section of the domain. Figure 5-21 shows some details of the computational grid in the wake of the diffuser wall region. In this calculation, a uniform velocity boundary condition is used at the inlet of the diffuser tubes. Boundary values for the k and e at the inlet are calculated based on 10% turbulence intensity and a length scale of 0.1 hydraulic diameter at the inlet. A uniform pressure boundary condition is used at the exit. To reduce the computational domain, a symmetry boundary condition is used in the region after the diffusers Plane area defined by 0.5375m<X <0.8375m, Y=0 and 0.5375m<X <0.8375 m , Y=0.0375m Computations are done using both the standard k-e model and Reynolds stress model. 84 Figure 5-19 : Computational geometry 85 0.4 h Figure 5-20: Computational grid. X Figure 5-21: Computational grid details in the wake region of the diffuser walls. To insure the grid independence of the numerical solution, three different grid arrangements used for numerical calculation for both turbulence models. Grid 1, 2 and 3 had 5300, 19200, and 66976 computational cells correspondingly. Figure 5-22 shows the 86 mean component of the velocity in the X direction for the region after the diffusers (along the line A B C in figure 5.19) for the three different grids. For simplicity only the results from the R S M models is shown. Figure 5.23 shows the turbulence kinetic energy at the same location as Figure 5-22 normalized with the 3/2 Ub 2 where Ub is the bulk velocity at the inlet to the converging section. Figure 5-23 shows similar results for the grid 2 and grid 3 indicating the grid independence of these two solutions. For simplicity only the results for grid 3 are presented in this section. Figure 5-22: Comparison of the mean velocity components in X direction in the region after the diffusers for different grid sizes. Velocities are nondimensionalized with the bulk velocity at the inlet to the converging section. X* =X/Xi where X t is the total length after the diffuser tubes . 8 7 Figure 5-23: Turbulence kinetic energy variations in the region after the diffusers for different grid sizes, k* is the turbulence kinetic energy nondimensionalized by the 3/2 (Bulk velocity)2 at the inlet to the converging section. Figure 5-24 shows the mean streamwise components of the velocity in the headbox measured along a vertical line passing through the point y=0.01875 m , x = 0.5675 m (3 cm after the diffusers exit). Computational values of the velocities obtained from both the k-e model and the R S M model show velocity profiles similar to the measured values which is different from the turbulent velocity profile that exists in the fully developed duct flow. However numerical values show a similar velocity profile for the top and bottom diffusers unlike the measured velocity profile which shows slightly different profiles for top and bottom diffusers. RMS values of the streamwise fluctuating component of the velocity at the same location as Figure 5-24 are shown in the Figure 5-25. 8 8 Figure 5-25 shows reasonable agreement between the measured and computational values obtained using the RSM model. Figure 5-26 shows the mean components of the velocity calculated numerically in the region after the diffusers at different locations along the center-line of the symmetry plane (line ABC in Figure 5-19). In Figure 5-26 U*, V* and W* are the mean components of the velocity nondimensionalized with the bulk velocity at the inlet to the converging section. There is reasonable agreement between the computed values of the mean velocities and their measured values. Figure 5-27 shows the variation of the turbulence kinetic energy along the same line as in the Figure 5-26. The turbulent kinetic energy predicted by both computational method in the regions immediately after the diffuser exits are very close to the measured values. However, the difference between the measured values and the computed values increases rapidly towards the exit. This shows that the turbulence models being used are not predicting accurately the actual turbulence in the contracting section. The RSM model is better than the k-e model, but even this more complex model is not adequate toward the exit of the contraction where accuracy is most needed if the fiber orientation leaving the slice is to be predicted correctly. The measured values of the turbulence kinetic energy in the converging section show a slight decrease which indicates that the production of the turbulence energy must be very small and in this way are somewhat similar to that of decaying turbulence flows where production is zero. Correct prediction of the turbulence kinetic energy in the headbox flow apparently requires the value of the production term in the turbulence kinetic energy equation to be close to zero. Although the standard k-e turbulence model performs well in simulations of homogenous decaying turbulence flows, such as grid turbulence where the production of turbulence is equal to zero, it fails to correctly predict the turbulence kinetic energy in the headbox flow. This indicates that an inaccurate estimation of the production of turbulence kinetic energy is being made by the k-e model in this case of significant plane rate of strain. This topic will be further investigated in the following chapter. 89 1E£ 0.9 0.8 0.7 0.6 ^ 0.5 0.4 0.3 0.2 F-0.1 -• v • • ke V RSM O Measurements 0.1 U (m/s) 0.2 Figure 5-24: Streamwise component of the velocity measured along the vertical line at the center of the diffusers at the location 3 cm after the diffuser exit plane. Z* is the height from the bottom plane divided by the total height at the location. 90 1 0.9 0.8 0.7 0.6 N 0.5 0.4 0.3 0.2 0.1 0 • • o Measu remen ts • R S M M o d e l I I I I I I I I I _L 0 0.01 0.02 X L 0.03 ff()4 0.05 u' rms (m/s) 0.06 0.07 0.08 Figure 5-25: RMS values of the fluctuating component of the velocity in the streamwise direction measured along the vertical line at the center of the diffusers at the location 3 cm after the diffuser exit plane. Z* is the height from the bottom plane divided by the total height at the location. 10r 8 * 6 JL 4 R S M U * A R S M W * o R S M V * o M e a s u r e m e n t s W * X M e a s u r e m e n t s U * * M e a s u r e m e n t s V * k-e U * • k-e W * ft k-e V * <3 <3 <I x •SI x x X X X X Figure 5-26: Comparison of the mean velocity components in the region after the diffusers. Velocities are nondimensionalized with the bulk velocity at the inlet to the converging section. X* =X/Xi where X| is the total length after the diffuser tubes . 9 2 6 -5 4 • k measurements < K-e 3 O RSM -<3 < 1 ' ' 1 1 < < < cP < < ] r P '0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X * Figure 5-27: Turbulence kinetic energy variations in the region after the diffusers. k* is the turbulence kinetic energy nondimensionalized by the measured k at the inlet to the converging section located at x =0.6185m(6 mm after the convergence starts). 5.3 Turbulence Kinetic Energy Increasing demands for higher paper quality including the need for control over the fiber angle in the final paper product, has put increased demands on the headbox design and ultimately requires the flow leaving the headbox to have a controllable level and degree of isotropy of turbulence. The average fiber angle leaving the headbox is governed by two competing factors, the mean flow shear tends to align the fibers in the machine direction and the turbulence eddies tend to randomize the fiber orientations [5]. Therefore, it is important for the computational model to simulate correctly the turbulence quantities. Poor computational predictions of the fiber orientation will almost certainly result from an unsuitable turbulence model. One of the incentives of this study 9 3 was to measure the detailed turbulence flow structure and use the measured values to select the most suitable turbulence model for the computations. In the following section, the turbulence quantities obtained from the different computational models are compared to the experimental results. From a practical point of view, at least two parameters are required to represent the turbulence flow characteristics: the turbulence kinetic energy, k, and the turbulence length scale. If the degree of isotropy is to be determined, equations for each individual component must be solved, involving further modeling constants and computational complexity. The computational geometry used in these simulations is the same as for the experimental model described in chapter 3. The computational geometry and boundary conditions are described in detail in chapter 5.2. The computational and experimental results is compared only for the region after the diffuser tubes since measurements were made only in this region. In the Figure 5-28 the rms data are normalized by the value of u' at the inlet (u' jn). Rms values of the fluctuating component of the velocity in the X direction, u' diminish towards the exit, while w' increases. Fluctuating components of the velocity in the Y direction v' remains approximately constant. Figure 5-28 shows that turbulence field becomes roughly isotropic a short distance inside the converging section. The turbulence field becomes very non-isotropic as the flow nears the exit. Turbulence models which are based on the assumption of isotropy, are not likely to give accurate estimates of the fluctuating components of the velocity or the overall kinetic energy, especially in the region close to the exit. Figure 5-28 shows the computed values of the u', v 'and w' normalized by the value of u' at the inlet (u'in) using the Reynolds stress model. There is a reasonable agreement between the numerical and experimental values. Different components of the Reynolds shear stresses normalized by the value of u'u'at the inlet (u'u'in) computed using the RSM are shown in Figure 5-29. The only significant component of the Reynolds shear stress is the u'w'. 94 3 3 2.5 h 2 F 1.5 1 0.5 o; -0.5 -1 4 • p • • • • • iT w' l{u'u')ln A v' w' /(u'u') l n O u' v' /(u'u'),,, • • • • • • • • • • 0 -i i i I 0.25 X* _l I l_ 0.75 _i I I L. Figure 5-29: Different components of the Reynolds shear stresses normalized by u'u' at inlet using the Reynolds stress model. 96 0.014 0.013 0.012 0.011 0.01 0.009 jo 0.008 "iE 0.007 0.006 0.005 0.004 ( 0.003 0.002 0.001 0 IVf pjaciirpmpntc A RSM Fluent 0> K-e with Production term=0 O Standard k-e ] • ^QD A X X 00 O o . A o A A O X X X X 1 1 1 1 1 ' 1 1 ' J 0.025 0.05 0.075 0.1 0.125 X(m) 0.15 0.175 0.2 0.225 Figure 5-30: Turbulence Kinetic energy at mid plane converging section. 1.1 1 0.9 0.8 0.7 CO J2 0.6 O • • A o o > Production of k standard k-e Epsilon standard k-e Production of K Measurments Production of k RSM Epsilon RSM • X(m) Fi term(epsilon) in the k equation for different models CR=10. In the previous chapter, it was shown that in general, the values of k from the standard k-e are much larger than the measured values. As has been noted already, it appears that, in the case of the flows with such a high longitudinal rate of strain, the production term in the k equation is not correctly modeled. Figure 5-30 shows the turbulence kinetic energy computed with different turbulence models compared with the experimental results. This figure shows a significant improvement when the k-e model in which the production term is arbitrary set to zero is used in the entire computational domain. To investigate further, the variation of the production and dissipation terms in the turbulence kinetic energy equation along the center-line of the converging section are shown in the Figure 5-31. It appears that the production term in the standard k-e model is highly over estimated compared to the similar term in the other turbulence models. The source of the error lies in the fact that, when using the standard k-e model, the production term is not calculated accurately. This can be shown by considering a simple 2D case where velocity accelerates in the X direction. 98 Figure 5-32: Geometry of a 2 dimensional converging section. For this flow the production term (equation 2.5) can be calculated as follows: P=-UiUjY.Sij where 5.. = fd U , d U j A + is the rate of stress tensor as defined in chapter 2. d Xj d X; Here: S,,=2 3U/ax, s22=2 av /ay For this case Si2=s2i=2(au/aY + av/ax)=o therefore p= - u'2 x 2 au/ax - v ' 2 x 2 av /ay p=2 (v ' 2 - M ' 2 )xau /ax (5.i) No assumption or simplification has been made to write equation (5.1) which is only a result of mathematical manipulation. This equation shows that the production terms in the k and e equations are proportional to the product of the dU/dX and the difference between w ' 2 andv ' 2 Since we do not calculate u' and v' explicitly in the k-e model, equation (5.1) cannot be used to calculate the production term in the turbulence kinetic energy equation. Therefore the k-e model approximates the terms v' and u' by relating them to the mean components of velocity. In the standard k-e model where the Boussinesq approximation is used, these terms are normally calculated as follows; 9 9 _ 1 uiuj--vijSij-~uiuidij so that, u'2 = -2/3 v, x 2 3U/3X - 2k2/3 5.2 v'2 = -2/3 v, x 2 av/aY- 2k2/3 = 2/3 v,x2 dU/dX - 2k2/3 These are clearly crude approximations, which are valid only when there are other stresses more important than the normal stresses (the usual case in turbulent shear flow for which the standard k-e turbulence model is designed). It should be noted that equation 5.2 makes the u negative when —>o. This is clearly impossible and shows the dx degree of inaccuracy of the standard k-e model in simulation of the flows with large positive au/ax . Therefore the only time that equation (5.2) can be used to approximate equation (5.1) is when the term dXJ/dX is small compared to other rates of strain such as dU/dY. In flows where the term dU/dX is the dominant term, the error between the actual value of the production term and its approximated value can become unacceptably large. However, if the flow is isotropic or very close to isotropic, then equation 5.1 implies that the production term is very small or zero. Since the measurements show that for most of the flow, except towards the exit, the turbulence in the headbox flow is close to isotropic, then the turbulence production should be very small or zero. Then for realism, one can set the production term to zero to produce better results than the standard k-e model. This has little impact on calculations of the mean flow quantities which are therefore estimated with reasonable accuracy. For the headbox flow where the term dU/dX is large and the shear rates of strain are small or zero, the error in the calculation of the k is unacceptable. To show this effect more quantitatively k values have been calculated for different contraction ratios (CR) and are plotted in the Figure 5-33 for the region after the diffuser tube along the line joining the center of the inlet and exit plane of the headbox region. The computational domain and the boundary conditions are the same as those described in section 5.2. In these calculations the total length in the streamwise direction 100 (x) was kept constant. Different contraction ratios are obtained by changing the angle between the top and the bottom plate. As expected, increased contraction ratios and therefore increasing 3U/8X values increase the difference between the k values obtained using the R S M and standard k-e model. CM E 0.02 x 0.06 -- • RSM CR=5 A RSM CR=10 -o RSM CR=20 k-e CR=5 0.04 — * k-e model CR=10 X k-e CR=20 0 0.05 0.1 0.15 X(m) X O o * X * 0.2 Figure 5-33: Turbulence kinetic energy k calculated for different contraction ratios. At best, the k-e model should be used to calculate turbulence kinetic energy k, only for small contraction ratios (CR < 5). 5.4 Turbulence Length Scale One of the important functions of the headbox is the prevention of formation and destruction of permanent floes in the pulp suspension. The study of the fiber flocculation and dispersion done by other researchers, especially the work done by Bennington et al. [60,61,62], suggests that in a fully turbulent regime of pulp suspension flow, fluid shear 101 entangles fibers to form a transient floe. This floe is often immediately dispersed. A permanent floe can only form when a transient floe does not disperse. In general, local turbulence affects the size and the strength of the transient floes. But the rate of decay of turbulence, and the floe velocity through this flow field determines whether the floe fiber entanglement remains transient or turns into a stable floe. The need for high fluid shear rate for dispersion of the fiber floes has led to widespread use of "high turbulence" headboxes. However, the degree to which turbulence can be increased in practice is limited by the need to maintain a stable jet emerging from the slice [60]. Grainy formation is another shortcoming associated with the high turbulence headboxes. Such formation is characterized by small floes having a high degree of contrast with the surrounding fibers in the paper. The rapid decay region, which usually follows the turbulence generating elements in the headbox, reduces the turbulence intensity. This reduction may allow reflocculation to set in before the paper is formed. In addition, the turbulence, if it is highly non-isotropic, can produce paper with properties which vary from one direction to another, due to the non-isotropic dispersion of fibers through turbulent action in the headbox. Therefore, in addition to modeling the turbulent kinetic energy, modeling the rate of decay and the degree of isotropy of turbulence are also necessary for evaluation of a headbox paper-making and deflocculation capability. Another turbulence characteristic, which has a great influence on the final paper quality, is the turbulence length scale [63]. In the past most researchers have used the assumption of isotropic decaying turbulence to estimate turbulence length scale and the rate of dissipation of the turbulence kinetic energy. Although the numerical simulation of the flow can be used to estimate these quantities, no study has been done in the past to show the range of accuracy of the different turbulence models in estimating these turbulence quantities in the headbox flow. In this chapter experimentally and numerically estimated turbulence length scales are compared. The rate of dissipation of the turbulence energy cannot be measured directly. However, it can be measured indirectly by measuring the turbulence dissipation length scale. In the past researchers have shown that for the case of the decay of turbulence behind a grid the turbulence integral length scale is very close to the dissipation length scale [64]. 102 Therefore the integral length scale can be used to estimate e in cases similar to decaying turbulence. The size of the turbulence integral length scale is dominated by the size of the larger, energy containing eddies associated with the turbulence. Direct measurement of the integral length scale requires simultaneous velocity measurements at two separate locations. However, it can be measured indirectly by measuring of the temporal time scale using single-point LDV method. In the presence of a strong uni-directional mean motion, Taylor's hypothesis can be used to relate the time scale to the length scale as follows: /= TxU c where Uc is the bulk velocity of the fluid and T and / are the integral time scale and turbulence length scale in the streamwise direction respectively. Temporal time scale are measured by integrating the autocorrelation function. The autocorrelation function at each measurement point is obtained using the FIND 403 software [53] which uses the lag slotted method to calculate the autocorrelation function allowing for selection of the lag slot time intervals. In the present study, only the streamwise component of the turbulent velocity is used to define the turbulent integral length scale. The computational length scales are calculated using 3 / = c „ x 4.4 £ where / is the turbulence dissipation length scale k and e are the turbulence kinetic energy and its rate of dissipation respectively. The computational domain and boundary conditions are the same as those described in chapter 5.2. In Figure 5-34, measured integral and calculated turbulent dissipation length scales are shown along the centerline of the mid plane (line ABC in Figure 5-19) for the region after the diffuser tubes. Due to insufficient LDV data rate, reliable experimental integral length scales could not be obtained in the region very close to the headbox exit. 103 As Figure 5-34 indicates, neither the trend nor the magnitude of the measured integral length scale is comparable to the calculated values the experimental values being between 3 times and one-half the calculated ones. As previously mentioned the size of the integral length is dominated by the size of the large, energy containing, turbulence eddies, while the dissipation length scale is affected by the turbulence kinetic energy and its rate of dissipation. Since the size of the large eddies is a function of the size of the model, one can conclude that the integral length scale is also strongly affected by the model size. In the headbox flow, while the characteristic length of the model, namely the height of the model, decreases the turbulence kinetic energy remains almost the same. That is why the integral length scale decreases as the height of the model decreases while the dissipation length scales remains almost the same. Therefore in headbox flow modeling the integral length scale of the turbulence should not be used to estimate the rate of dissipation of turbulence kinetic energy. Figure 5-34: Variation of the length scale along the centerline in the mid plane near the entrance. L * is the length scale divided by the width(=0.0375m) of the diffusers at the exit. X* is the distance in the x direction after the diffuser tubes divided by the total length of the headbox after the diffuser tubes 104 5.5 Free surface Jet In this section, the computation domain is extended to include the free surface jet beyond the headbox. The VOF model described in chapter 2 was used for the numerical simulation of the flow through contracting section of a headbox and the jet after the headbox. Figure 5-35 shows the computational domain and its respective dimensions. Boundary conditions are as follows At inlet to the converging section(X=0), for plane area : X=0 , 0<Y<0.075 and 0.1<Z< 0.175 Uniform velocity boundary condition U=0.2 m/s =Constant k= 0.0035 m2/s2 e=0.011m2/s3 Water phase volume fraction aw=l For Plane area: X=0 , 0<Y<0.075 and 0<Z< 0.1 and 0.175<Z<0.275 Pressure Boundary condition P = P 0 =Constant where P 0 is the reference pressure =101 kPa Air phase volume fraction ota^ l At the exit for the plane area: X=0.325, 0<Y<0.075 and 0<Z<0.275 Gradient Boundary condition dU/dx = 0 At the top for the plane area: 0<X<0.325, 0<Y<0.075 and Z=0.275 Pressure Boundary condition P=P0 At the bottom for the plane area: 0<X<0.325, 0<Y<0.075 and Z=0 Pressure Boundary condition P=P0 On the south side for the plane area: 0<X<0.325, Y=0 and 0<Z<0.275 Symmetry boundary condition 3V/3Y=0 105 On the south side for the plane area: 0<X<0.325, Y=0.075 and 0<Z<0.275 Symmetry boundary condition dV/dY=0 A typical computational grid is shown in Figure 5-36. Grids are concentrated near the air-water interface to reduce the numerical diffusion and to increase the resolution for tracing the interface. The SIMPLE (Semi Implicit Method for Pressure-Linked Equations) was selected for the solution algorithm. A second order up wind was chosen for the solution of all the transport equations. Both standard k-e and Reynolds stress model was used for modeling the turbulence quantities. Figure 5-37 shows the variation of the mean streamwise component of the velocity in the mid plane along the line connecting the center of the inlet plane to the center of the exit plane (line AB on the Figure 5-35). After the exit from the headbox numerical data on Figure 5-37 are taken along the center-line of the jet. In the Figure 5-37 the U* is the mean streamwise component of the velocity divided by the bulk velocity of the water phase at the inlet to the headbox section and X* is the distance along the X-axis divided by the total length of the converging section. Figure 5-37 shows that the U* values do not decrease immediately after the exit plane indicating the existence of a short contraction after the exit. Similar to the internal headbox flow simulations there is a reasonable agreement between the experimental and numerical results (using both standard k-e and Reynolds stress model) as far as mean flow calculations are concerned. Figure 5-38 shows the root mean square values (rms) of the velocity fluctuations normalized by the rms value of the velocity fluctuation in the streamwise direction (u' jn) at the inlet at the same positions as Figure 5-37. Overall the variations of individual normal Reynolds stresses for the region inside the headbox are similar to those shown in Figure 4-7. The streamwise component u 1 decreases while the Z direction component w 1 increases and Y direction (cross machine direction) component v 1 stays approximately the same resulting in a non-isotropic turbulence field especially toward the exit. However this trend changes after the exit from the headbox. The w1 decreases rapidly in value and the turbulence field moves toward an isotropic condition. Unfortunately there are no experimental data available close to the exit for the region inside the headbox but in the 106 regions where the experimental results are available, there is a reasonable agreement between experimental and computational results. Figure 5-39 shows the turbulence kinetic energy variation at the same location as Figure 5-37. In Figure 5-39 X* is the distance in the X direction normalized by the length of the converging section and the K* is the turbulence kinetic energy k normalized by its value at the inlet. Since there were no experimental values of w' available the experimental values of the k are calculated assuming that w' is equal to the average of u' and v'. Examination of Figure 5-38 shows some inaccuracy of this assumption, but this result is included in the figure to give an estimate of the experimental k. Numerical (using both k-e and Reynolds stress model) and experimental results show a rapid reduction in turbulence kinetic energy in the jet after the headbox exit. Pictures taken of the jet in the experimental set up (Appendix B) show ripples on the surface of the jet. These ripples are absent in the numerical simulations partly due to the fact that these simulations are based on the steady state assumption. The numerical simulations of this chapter show that after the initial decay of the turbulence there is short region close to the exit inside the headbox where actually the turbulence kinetic energy increases. However the decay of the turbulence energy is present in the jet exiting from the headbox. The results from the k-e model in the jet have the same trend as the results from RSM model as shown in Figure 5-39. However, the results from k-e model stay at unrealistically high values through the jet. This is due to the error that was introduced in the k-e results from computations inside the converging section of the headbox. Although the k-e model usually performs reasonably well in predicting free surface jets, it should be used cautiously in the calculation of the headbox flow, since the error that originates from the calculation inside the converging section will continue to affect the computational results in the jet as well. However, if the k-e model is used in the free surface jet, the calculation should be started with the boundary values at the start of the free surface extracted from computation of the flow inside the headbox using RSM or the k-e model with the generation term neglected. Figure 5-40 shows the variation of the Reynolds shear stresses at the same locations as Figure 5-37. In Figure 5-40 Reynolds shear stresses are normalized with the 2/3 xk at the inlet to the headbox section. The only significant Reynolds shear stress is the u'w 1 107 component which is negative at the beginning but rises toward the exit of the headbox followed by a rapid decrease in the free surface jet. Figure 5-35: Free surface jet calculation domain. 108 Figure 5-36: Numerical Grid used for free surface calculation. p' 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 X* Figure 5-37: Mean component of the velocity in the streamwise direction nondimensionalized with the inlet bulk velocity. X* is the distance along the X axis normalized with the total length of the converging section 109 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 : • RSM model u'rms : A RSM model v'rms * RSM model w'rms : Measurements u'rms ~ 0 Measurements v'rms I 0 0.25 0.5 0.75 Figure 5-38: Rms. values of the fluctuating component of velocities u 1 , v 1 ,w normalized by u 1 inlet. 2.5 2 . 1.5 1 0.5 0 - • RSM model - ke model - o Experiment V V V v V V v v v v v v v v v v v v W • 0 0.5 Figure 5-39: Turbulence kinetic energy along the centerline of the mid plane inside the headbox and along the center of the free surface jet outside the headbox. X* =1 is the exit -plane of the headbox. K* is the turbulence kinetic energy normalized with k value at the inlet. 0.4 0.35 0.3 3 0.25 0.2 0.15 • uw* A vw* o uv* • • • • • • • Figure 5-40: Variation of Reynolds shear stresses in the free surface jet. Reynolds shear stresses are normalized with 2/3 k at the inlet. X* is the distance along the X axis divided by the total length of the headbox. i l l 6 - Summary and conclusion In this thesis, the flow through a paper machine headbox has been investigated both experimentally and numerically. Particular attention was given to turbulence characteristics and turbulence generated secondary flows. Experimental measurements of the mean and fluctuating components of the velocity were carried out in a small scale down model of a typical headbox. These measurements provide a reliable test case for validation of numerical methods to be used for modeling flow in a headbox. Different numerical models were used to simulate the complex turbulent flow in a paper machine headbox. Numerical results were compared to the experimental results and the accuracy of each turbulence model in simulating the mean flow and turbulence quantities was examined. These numerical results are complementary to experimental measurements by providing detailed information of the flow where experimental measurements are not possible. The major findings of this thesis are summarized as follows. 6.1 Summary 6.1.1 Experimental measurements LDV measurements of the mean and fluctuating components of the velocity at different location inside a paper machine headbox model were done. Mean flow measurements show that flow inside the headbox is almost two-dimensional except for small regions close to the walls. These measurements showed that although the footprint of diffuser tubes vanishes after a short distance into the converging section, there is slight variation of the velocity profile in the paper thickness direction even close to the exit. Measurements of the three fluctuating components of velocity showed that they have different values at the start of the converging section with the streamwise component having the highest value. Streamwise rms values of the velocity fluctuation u 1 decrease toward the exit of the converging section while the w"(paper thickness direction) increases and the v' (cross machine direction) stays approximately the same. Turbulence 112 anisotropy decreases initially but starts increasing a short distance inside the headbox particularly near the exit. Turbulence integral length scale was also measured experimentally showing that the turbulence integral length scale decrease rapidly towards the exit, reflecting the decaying size of the headbox channel. Measurements in the jet emanating from the headbox show small contraction of the jet after exiting from the headbox. These experimental results also show a rapid reduction in turbulence kinetic energy in the jet after the headbox exit. 6.1.2 Numerical simulation of the converging section 6.1.2.1 Boundary conditions Different boundary conditions greatly influence the accuracy of the numerical results. Numerical simulation of the flow in the converging section of the headbox was carried out with the application of different exit boundary conditions, namely pressure and modified gradient boundary conditions at the exit. For uniform flow entering the converging section these two boundary conditions produce similar results. However, the results from these two boundary conditions at the exit differ from each other significantly if the flow at the inlet to the converging section is not uniform. In the case of non-uniform inflow it appears that the pressure boundary condition produces more accurate results. Flow in symmetric and non-symmetric converging sections were compared. For the case of uniform flow at the inlet, the flow pattern and the magnitude of the secondary flows at the exit are the same for both symmetric and non-symmetric converging section. However for the non-uniform inlet condition these two geometries produce different flow patterns. A cork screw motion is produced in the non-symmetric converging section. This cork screw motion results in strong cross flow velocities at the exit. In general symmetric converging sections are more efficient in reducing the effects of the non-uniform inflow than non-symmetric converging sections. 113 6.1.2.2 Secondary flow investigation The secondary flows, which are in the cross machine direction, are very small compared to the machine direction velocity but they can become important as soon as the jet from the headbox touches the moving fiber mat. This is due to the fact that the fiber mat is moving approximately at the same speed as the jet therefore the velocity of the stock relative to the fiber mate in the MD becomes nearly zero and the CD velocity becomes the only significant component of velocity which can effect the fiber orientation angle. Numerical simulations of the flow in the headbox were carried out using different turbulence models to investigate the turbulence generated secondary flows. The turbulence models used were the standard k-e model, quadratic and cubic algebraic turbulence models and the Reynolds stress model. The results showed that turbulence generated secondary flows are small in magnitude and are confined to a very small region close to the wall. In most practical cases the presence of the non-uniformity that may originate from any unbalanced flow in the distributor will be much more important, and these turbulence generated flows can be ignored without significant loss of accuracy. The results from both quadratic and cubic algebraic turbulence models showed a significant improvement over the results from standard k-e in modeling the turbulence generated secondary mean flows. However, they did not show much improvement over the standard k-e model in estimating the turbulence kinetic energy in the converging section of the headbox. 6.1.2.3 Full headbox model The numerical modeling of the full headbox was carried out using the standard k-e model and the Reynolds stress model and the results were compared to experimental measurements. Standard k-e model and Reynolds stress model predicted the mean measured velocity components reasonably well. Comparison of measured rms values of the velocity fluctuations with the values from the Reynolds stress model show good agreement between the numerical and experimental results. Other turbulence models were less successful. 114 6.1.2.4 Turbulence kinetic energy In the results of a simulation of the headbox flow using the standard k-e model, the turbulence kinetic energy values were much higher than those predicted from the Reynolds stress model or the experimental values. The difference between the standard k-e and Reynolds stress model values for turbulence kinetic energy k increases as the contraction ratio is increased. In general, the standard k-e model values are close to the values obtained from the Reynolds stress model only for contraction ratios less than 5. A significant improvement in the results from the k-e model can be achieved by artificially setting the production of the turbulence energy to zero. 6.1.2.5 Turbulence Length scale Values of the turbulence integral length scale were measured using the LDV "Slot" technique. Comparison of the computed dissipation length scale with measured integral length scale showed disagreements between numerical and experimental data both in trend and the values. Therefore, measured integral length scale should not be use to predict the rate of dissipation of the turbulence kinetic energy e. 6.1.2.6 Free surface jet The volume of fluid method was used to model the three dimensional flow in the headbox and the free surface jet exiting from the headbox. The standard k-e and the Reynolds stress models were used for turbulence modeling. Comparison of the mean velocity components showed good agreement between the numerical results (both standard k-e and Reynolds stress model) and the measured values. The streamwise velocity component profile shows a small contraction of the jet after the headbox. There is a good agreement between the values obtained by using the Reynolds stress model and the experimental measurements for the fluctuating components of the velocity. After the headbox, in the free surface jet the w1 decreases in value and the turbulence field moves 115 toward an isotropic condition. The only significant Reynolds shear stress is the u' w' component, the other two shear stresses being very close to zero. Measurements and numerical values from the Reynolds stress model show that after an initial decaying of the turbulence near the entrance of the contraction, there is a short region close to the exit inside the headbox where the turbulence kinetic energy actually increases. Decaying of the turbulence energy continues in the jet exiting from the headbox. The results from the k-e model in the jet have the same trend as the results from RSM model. However, the results from k-e model stay at unrealistically high values through the jet. This is due to the error introduced in the k-e results from computations inside the converging section of the headbox. Although the k-e model usually performs well in predicting free surface jets, it must be used cautiously in the calculation of the headbox flow, since the error that originates from calculation inside the converging section will substantially affect the computational results in the jet. In general the results of these study can be used to improve the design of the headboxes to provide a greater control over the average fiber angles in the jet emerging from the headbox. The fiber angle in the jet is determined by two competing phenomenon in the headbox: mean flow normal rate of strain, which tends to align the fibers in the machine direction and the turbulence fluctuation which tends to randomize the fiber orientations. Usually the mean flow is determined by the need for it to have a value very close to the machine speed, which in turn is determined by the paper production rate. However the turbulence field can be varied to control the fiber angle. By knowing the turbulence field the length of the converging section can be varied so that the turbulence kinetic energy is, for example, at its lowest value close to the exit which will line up the fibers in the machine direction. If the length of the converging section is selected so that the kinetic energy is at a higher value close to the exit, the fiber angles will be more randomly oriented. For existing headboxes with fixed length, external forces (such as electromagnetic or ultrasound) can perhaps be used to change the initial values of the turbulence kinetic energy which in turn would change the fiber angle orientation in the jet. 116 The numerical modeling guidelines provided in this thesis could be used to model the flow in any headbox geometry and so could be used to estimate the turbulence quantities and thereby to predict the fiber angle in the jet from the slice. Numerical models of the entire headbox are possible and the present work provides further developments of the accuracy of existing methods. 117 6.2 Conclusion Mathematical modeling of headboxes can contribute very significantly to the understanding, design and operation of the headboxes at reduced cost. The model is very complex and has to been done with a thorough understanding of the physics. Some of the conclusions and recommendations for developing mathematical models are listed below. • Turbulence generated secondary flows in the headbox are not as significant as reported earlier by previous authors and their effects on the main flow are limited to a small region close to the side walls. • The main flow is predicted with sufficient accuracy by turbulence models used in the present study namely the standard k-e model, quadratic and cubic algebraic turbulence models and Reynolds stress model. • Local tube footprint (flow profile at the exit of each diffuser tube and the wake generated by the wall separating the diffuser tubes) tends to be small in the mean flow and diminishes rapidly as the flow enters the converging section. Non-uniformities originated from the header tend to be preserved in the convergent section. • The turbulence predictions for high contraction ratios (CR > 5) using turbulence models based on Boussinesq approximation (such as standard k-e model) are not accurate. These predictions are somewhat better for CR < 5 and in these cases, turbulence models such as the standard k-e, can be used for turbulence predictions with reasonable accuracy. • The lack of performance of k-e is due to over estimation of the generation term in the turbulence kinetic energy equation. Therefore ignoring the generation term (arbitrarily setting the generation term equal to zero) actually improves the predictions of the turbulence quantities. • Turbulence kinetic energy is relatively constant through the slice and increases somewhat towards the exit. However, relative turbulence intensity is much lower closer to the exit of the converging section. Turbulence also becomes significantly anisotropic with the Z (paper thickness direction) 118 component being much larger than the other two components, in agreement with well understood models of the effect of distortion on turbulence[64]. This trend will affect the distribution of fibers at the headbox exit and should be accounted for when modeling fiber distribution. In numerical modeling of headbox flow, the pressure boundary conditions should be used at the exit of the slice, particularly when the inflow of the contracting section is not uniform. The gradient condition appears to be inaccurate in the case of the non-uniform inflow. The free surface jet after the exit of the slice can be computed with k-e. However, if the standard k-e used for the calculation inside the slice the results are erroneous as the computed turbulence is too high at the exit of the slice and this error is carried into the free jet. 119 7 - Recommendation for future work Further numerical and experimental works can be done to extend the knowledge of the headbox and its modeling. The primary area of the experimental work should be the study of the free surface jet after the headbox. High speed photography and automatic digitization and analysis of the surface ripples would give an insight into the origin and development mechanism of the surface ripples. Statistical analysis of the free surface jet with the different flow arrangements inside the headbox could be beneficial in determining if the diffuser tube size and arrangement have any effect on the size and the location of the surface ripples. In addition, an experimental study of stratified headboxes, where different fluids layers are separated with thin flexible membranes, would give a better understanding of the turbulence development and mixing inside these headboxes. The effects of the lip profiling of the headbox which were not considered in this study, should be examined both numerically and experimentally to understand the effect of this geometry change on the turbulence field and mean flow. Numerical simulation of the headbox flow using the complex Large Eddy Simulation (LES) model can further expand the work of this thesis. LES can provide a better understanding of the turbulence field inside the headbox and can provide more detailed information in the areas where experimental measurements proved to be difficult (near the wall and close to the exit). The volume of the fluid method or another similar numerical method capable of simulating the free surface should be used in transient mode to capture the surface ripples. This will increase the usefulness of CFD as a tool for study of the headbox design. 120 Nomenclature , Ci, C 2 , CD, CE Turbulence Model Constants C 3 , C 4 , C 5 , C6, C 7 Turbulence Model Constants CD Cross machine direction dj Convection term in the Reynolds stress equation CR Contraction ratio Dh Hydraulic diameter DjjL Molecular Diffusion term in the Reynolds stress equation DyT Turbulent diffusion term in the Reynolds stress equation E Wall function constant Production by system rotation in the Reynolds stress equation Gy Buoyancy Production term in the Reynolds stress equation L Turbulence length scale N Number of points in the data P Production rate of turbulence kinetic energy Py Stress Production term in the Reynolds stress equation Re Reynolds number T Turbulence intensity Sy Rate of stress tensor U,V,W Mean velocity in the X,Y and Z direction, respectively Up Velocity component parallel to the wall at the node p de Diameter of beam at probe(LDV) 121 df Fringe spacing d m Diameter of the measurement volume fs Frequency shift lm length of measurement volume i T k Turbulent kinetic energy (= 1/2 M -M . ) p Modified pressure term t Time y + Dimensionless normal distance from the wall yp Normal distance from the wall u', v', w' Turbulent (r.m.s) normal stresses uT Friction velocity Greek symbols a Volume fraction coefficient a, p\ y Alignment angles in x-,y- and Z direction respectively e Rate of dissipation of turbulent kinetic energy per unit mass 6ij Dissipation term in the Reynolds stress transport equation § General transport quantity <|)jj Pressure strain K Half angle of the beam pair (LDV) A, Wavelength \i Molecular viscosity p Density 122 0\, 0"e Turbulent model constant Subscripts Reynolds stress tensor (= puu • ) Wall shear stress Kinematic viscosity 8 Refers to time for eight fringe crossing for LDV Fringe spacing i n Values at the inlet Superscripts Turbulence quantity Nondimensionalized value Bibliography 1. Newcombe, D. "The headbox : First link in the quality chain", PIMA : 31-34, March (1990). 2. Kallmes, O. and Thorp, B. , " Back to basic in headbox design and application", In TAPPI Papermakers conference: 221 -226, (1984). 3. Smook, G. A., Handbook of Pulp and Paper Technology, TAPPI, (1982). 4. Britt, k. W., Handbook of pulp and paper technology. Van nostrand reinhold company, (1993). 5. Zhang, X. ,"Fiber orientation in Headbox.", Master Thesis, University of British Columbia, (2001). 6. Hua, L., "Computational Modeling of a Manifold Type Flowspreader", M.A.Sc. Thesis, University of British Columbia, 1998. 7. Mardon, J., Hauptmann . E.G., Monahan. R. E. and Brown. E. S., "The Extant State of the Manifold Problem", Pulp &Paper Magazine of Canada, 72(11):76-8, (1971). 8. Mardon, J., Monahan. R. E. and Brown. E. S., "Experimental Investigation of Manifolds Systems for Paper Machine Headboxes", Pulp & Paper Magazine of Canada: 459-474,( Sept 1965). 9. A.D. Trufitt. Design Aspects of Manifold Type Flowspreaders. In Pulp and Paper Technology Series, TAPPI,(1975). 10. Syrala, S., Saarenrinne, P. and Karvinen, R.," Fluid Dynamics of a Tapered Manifold Flow Spreader", In TAPPI Engineering Conference: 223-23, Chicago, Illinois, (1988). 11. Jones G. L. and Ginnow. R. J.," Modeling Headbox Performance with Computational Fluid Mechanics.", In TAPPI,Engineering Conference: 15-20, Chicago, Illinois, (1988). 12. Shimizu. T. and Wada, K., " Computer Simulation of Measurement of Flow in a Headbox.", In Proceedings of the,Pan Pacific Pulp and Paper Technology Conference: 157-165, (1992). 13. Hamalainen, J.," Mathematical Modelling and Simulation of Fluid Flows in the Headbox of Paper machines." Ph.D. Thesis, University of Jyvaskyla, (1993). 124 14. Tarvainen, P., Makinen, R. and Hamalainen, J.,"Shape Optimization for Laminar and Turbulent Flows with Applications to Geometry Design of Paper Machines.", Proceedings of the 10th International Conference on Finite Elements in Fluids, University of Arizona: 536-545, (1998). 15. Lee, J. J. and Pantaleo , S. B.," Headbox Flow Analysis.", In 84"" Annual Meeting, Technical Section CPPA: B339-B344, (1998). 16. He, P., Bibeau, E. L., Hua, L , Salcudean. M. and I. Gartshore, " Fluid Dynamics of the Flow Distribution in a Headbox", CPPI Conference, Montreal, (January 1998). 17. Hua, L, Bibeau, E. L., He, P., Salcudean. M. and I. Gartshore, " Flow Distribution in a Hydraulic Headbox." TAPPI Conference, Anaheim, (Sept 1999). 18. Hua, L, Bibeau, E. L., He, P., Salcudean. M. and I. Gartshore," Three Dimensional Flow Distribution in a Headbox Manifold", CANCAM99 conference, Hamilton, Ontario, (April 1999). 19. Hua, L., He, P., Salcudean. M., Gartshore. I. and Bibeau, E. L . , " Turbulent Flow in Hydraulic Headbox", TAPPI Conference, Vancouver, (April 2000). 20. Aidun, C. K. and Kovacs, A., "Hydrodynamics of the Forming Section: The Origin of Non-uniform Fiber Orientation", TAPPI 1994 Engineering Conference, Atlanta GA, (1994). 21. Aidun, C , "Hydrodynamics of Streaks on the Forming Table." TAPPI Journal, 80(8): 155-162, (1997). 22. Bandhakavi, V. S. and Aidun, C. K., "Analysis of Turbulent Flow in the Converging Zone of a Headbox", 1999 TAPPI Engineering/Process and Product Quality Conference & Trade Fair: 1135 (1999). 23. Parsheh, M., Dahlkid, A. /'Numerical Modeling of mixing in stratified headbox jet", 1997 Engineering and Papermakers conference, Nashville, (1997). 24. Parsheh, M. /'Aspects of boundary layers, turbulence and mixing in plane contraction with elastic vanes, investigating potential application to paper manufacturing.", Ph.D. Thesis, Royal Institute of technology, Stockholm, Sweden. (1999). 125 25. Zhang, X., Gartshore, I. and Salcudean, M.," Fiber orientation in a converging headbox flow field: measurements and numerical prediction", Submitted to Journal of pulp and paper Science, (2001). 26. Ullmar, M. and Norman, B.," Observation of fiber orientation in a headbox nozzle at low consistency ", Tappi Conference, Nashville, October (1997). 27. Ullmar, M., " On fiberorientation mechanisms in a headbox nozzle.", Master's thesis, Royal Institute of Technology, (1998). 28. Ferrier, A. C. and Aidun, C , "High-speed digital imaging of paper forming: A method for quantitative evaluation of paper forming hydrodynamics.", TAPPI Conference, Nashville, (October 1997). 29. Aidun, C , "Quantitative evaluation of the forming jet delivered from four different hydraulic headboxes using high-speed digital imaging.", TAPPI Conference, Nashville, (October 1997). 30. Soderberg, D. /'Experimental and theoretical studies of plane liquid jets.", Ph.D. thesis, Royal Institute of technology, (1997). 31. Soderberg, D. and Alfredsson, P.H., "Experiments concerning the creation of streaky structures inside a plane water jet.", TAPPI Conference, Nashville, (October 1997). 32. Nowak, P., "GEO-MGFD: A system of computer Programs for Steady-State Flow of Variable Density Fluids in the Complex Three-Dimensional Regions using Multi-grid and Multisegment Techniques", Technical Report, Dept. of Mech. Engineering university of British Columbia, Vancouver, Canada, (1992). 33. He, P. and Salcudean, M., "A Numerical Method for 3D Viscous Incompressible Flows Using Non-Orthogonal Grids." International Journal of Numerical Methods, in Fluids 18: (1994). 34. He, P.," Numerical prediction of film cooling of turbine blades using multi-block curvilinear grids.", Ph.D. thesis, University of British Columbia, Vancouver, Canada, (1995). 126 35. Boussinesq, J., " Theorie de l'ecoulement tourbillant", memoires presentes par divers savants sciences mathematique et physiques, Academie des sciences, Paris 23, (1877). 36. Launder, B. E. and Spalding, D. B., "The Numerical Computation of Turbulent Flows", Computer Methods in Applied Mechanics 3: 269 (1974). 37. Launder, B. and Spalding, D., Mathematical models of turbulence. London: Academic press, (1972). 38. Launder, B. and Spalding, D., " A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows", Int. J. Heat and Mass Transfer, (15): 1787-1806,(1972). 39. Speziale, C. G., "On nonlinear k-1 and k-e models of turbulence", Journal of Fluid mechanics, (178): 459-475, (1987). 40. Speziale, C. G , Gatski, T. B. , "On explicit algebraic stress models for complex turbulent flows.", Journal of fluid mechanics,( 254): pages 59-78, (1993). 41. Launder, B. E., Craft, T. J. and Suga, K.," Extending the applicability of eddy viscosity models through the use of deformation invariants and non-linear elements.", Proceedings of 5th international symposium on refined flow modeling and turbulence measurements: 125-132, Presses ponts et Chaussees, Paris, (1993). 42. Hoagland, L. C , " Fully developed turbulent flow in straight rectangular ducts: secondary flow , Its cause and effect on primary flow", Ph.D. Thesis, MIT, (1960). 43. Gessner, F. B.," Turbulence and mean Flow Characteristics of Fully Developed Flow in the Rectangular Channels." ,Ph.D. Thesis, Purdue university , (1964). 44. Gessner, F. B., Po, J. K, and Emery, A. F.," Measurements of Developing Turbulent Flow in a Square Duct" 1st international symposium on Turbulent Shear Flow. Perm. State University (1997) 45. Launder, B.E., and Ying, W. M . , " Secondary Flows in Duct of Square cross-section" Journal of Fluid Mechanics, (54): (1972) 46. Brundntt, E, and Baines, W. D., "Production and Diffusion of Vorticity in Duct Flow", Journal of Fluid Mechanics, (19): (1964) 47. Launder, B.E., and Ying, W. M. "Prediction of Flow and Heat Transfer in Duct of Square Cross-Section." Heat and Fluid Flow, 3, (2): (1973). 127 48. Nakayama, A., Chow, W.L., and Sharma, D., " Calculation of Fully Developed Turbulent Flows in Duct of Arbitrary Cross-section.", Journal of Fluid Mechanics, (128): (1983). 49. Demuren, A. O., and Rodi, W.," Calculation of Turbulence Driven Secondary Motion in Non-Circular Duct." Journal of Fluid Mechanics, (140): 189-222, (1984). 50. Myong, H.K., Kobayashi, T., "Prediction of Three Dimensional Developing Turbulent Flow in a Square Duct with an Anisotropic Low Reynolds number k-e Model.", Journal of Fluid Engineering ,( 113): 608-615, (1991). 51. White, F. M., Fluid Mechanics fourth edition. Mc Graw Hill. (1998). 52. Durst, F., Melling, A. and Whitelaw, J., Principles and Practice of Laser Doppler Anemometry (2 ed.). Academic press, (1981). 53. TSI Inc., Flow Information Display (FIND) Software Version 4.0 Instruction manual, TSI Inc., (1993). 54. Nakao, S., Terao, Y. and Hirata, K. ," New Method for eliminating the statistical bias in highly turbulent flow measurements", AIAA Journal 25(3):443-447, (1987). 55. George, W. and Lumley, J.,"The laser Doppler velocimeter and its application to the measurements of turbulence.", Journal of Fluid Mechanics 60: 321-362, (1973). 56. Edwards, R. /'Report of the special panel on statistical particle bias problems in laser Anemometry", ASME Journal of Fluids Engineering 109: 89-93, (1987). 57. Fuchs, W., Nobach, H. and Tropea, C. ," The simulation of LDA Data and its Use to Investigate the Accuracy of Statistical Estimators", Technical paper submitted to AIAA Journal: 1-21,(1993). 58. Kline, S. J.," The purposes of uncertainty analysis", Journal of Fluids Engineering, 107:153-160,(1985). 59. Findly, M. " Experimental and Computational Investigation of Inclined Jets in a Crossflow", Ph.D. thesis, University of British Columbia, (1998). 60. Bennington, C . , Kerekes, R. J.," Power requirements for pulp suspension fluidization", TAPPI Journal, 78(2) :253-258, (Feb. 1996). 61. Bennington, C . , Kerekes, R. J., Grace, J.R.," Motion of pulp fiber suspensions in rotary devices", Canadian Journal of Chemical Engineering 69(1): 251-258 (Feb. 1991.) 128 62. Bennington, C . , Kerekes, R. J., Grace, J. R.," The Yield stress of Fiber Suspensions", Canadian Journal of Chemical Engineering 68(5): 748-757 (Feb. 1990) 63. Olson, J. A. and Kerekes, R. J., " The motion of fibers in turbulent flow", Journal of Fluid Mechanics, 377: 47-64, (1998). 64. Townsend, A. A . , The Structure of Turbulent Shear Flow. 2n d edition, Cambridge University press, (1976). 129 Appendix A Angular Alignment Measurements In order to determine the angles of alignment of the probes with respect to the model coordinate system 10 measurements of each angle at various locations in the experimental setup were made and the standard deviation of the angles was calculated. The traverse mechanism and the probe angle were adjusted till the standard deviation in each angle was below 0.05 degree. The angles about the x, y, and z-axes are referred to as a, (3 and y respectively, as shown in figure A. 1. The measurement of the angle a will be described. The measurement of the angle p is similar. A different method was needed to measure y, which will be presented. A similar approach is used for measuring these angles when the probe was used to measure the z component of the velocity. To determine a a machinist's square was placed on the bottom plate of the headbox model. Using the unshifted component of the green beam pair initially, the y-traverse was adjusted until the beam intersected a mark on the arm of the square, which was perpendicular to the floor. The position of the y-traverse at this location was recorded as yi. The y-traverse was then moved which resulted in the point of beam intersection moving along the square until another reference mark was reached. The separation of the two marks, Az, then was known and the new position of the y-traverse was recorded as y2 . The value of a for the green beam was calculated from: a = arctanl To determine y the pinhole block was used [59]. The shifted green (or blue) beam was positioned over the pinhole and the positions of the traverse mechanism were recorded (xi, yl). The traverse was then adjusted in the x and y directions until the unshifted green (or blue) component was centered over the pinhole where X2 and were recorded. The angle was calculated from: 130 y = arctanl I v 2 - v , ; where the yp is the angle obtained using the pinhole method. The angle Yp is not the true angle of alignment as there may be an error from the traverse itself if the traverse directions are not perfectly parallel to the corresponding tunnel directions. To measure this angle, the y-traverse (in the case of the green beam pair) was moved over a known distance (from the traverse), Ay, along a line parallel to the y-axis in the wind tunnel. The x-traverse was adjusted at the end of the traverse, if necessary, to bring the beam back in line with the y-axis. The locations at the beginning and end (xi and X2) of the traverse were recorded and the traverse angle was calculated from: / y, = arctan -Ay where the Yt is the traverse angle. The angle of alignment was then: Y = YP-Yt A similar procedure was used to calculate and adjust the angles when the probe was used to measure the z direction component of the velocity. 131 z 7 ° 1 P x Figure A.1: Alignment angles 132 Appendix B Free surface ripples Pictures of the free surface jet were taken during the course of this study to show the instability of the jet exiting the headbox. Pictures show the existence of ripples and non uniformities on the jet surface that do not appear in the numerical simulation. These pictures were taken using a Sony DCR-TRV320 digital video camera mounted in front of the headbox model pointing in the - X direction at the same height as the exit of the headbox. The lighting was provided from the top. The Sony camera has a shutter speed of 1/4 to 1/4000 second. Figure C.l shows the schematic of the camera location. Observation of the jet showed two different kinds of surface waves: stationary waves and dynamic or moving waves. Stationary waves in the form of a "V" originated from surface irregularities at the edge of the headbox walls. The population and the depth of these waves were reduced greatly by sanding and smoothing the edges of the headbox model. The jet exiting from the headbox had a smooth surface for a short distance after the exiting from the headbox where only stationary waves were apparent. This smooth region was followed by the appearance of dynamic waves, which grew in amplitude with no apparent fixed location. Light Source Camera Figure B-l: Camera Location 133 Stationary Figure B.2: Stationary and dynamic waves: Picture B.2 shows example of the stationary and the dynamic waves already described. The pictures taken at a higher shutter speed than the one used in figure B.2 show some type of cross-wave appearing in the surface of the jet. The appearance of this type of waves was also reported by Soderberg [30]. Picture B.3 shows a typical picture of these cross-waves. 134 Figure B.3: Cross-waves in the free surface jet 135 Lines corresponding to the location of the tube walls upstream Dynamic waves Figure B.4: Location of the Dynamic waves. From visual inspection of the pictures it appears that neither the amplitude nor the location of the dynamic waves are related to the width of the tubes; these visual inspections are not conclusive and this observation requires further investigation. The relationship between wave crest location and tube position requires an experimental set up where many pictures can be taken and then digitized and analyzed automatically to find the density function distribution of the wave crests. 136
Thesis/Dissertation
eng
Mechanical Engineering
Vancouver : University of British Columbia Library
For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use.
Graduate
Expermental and mathematical modeling of flow in headboxes
Text