HYDRODYNAMICS OF THE HYDROCYCLONE FLOW FIELD: EFFECTS OF TURBULENCE MODELLING, GEOMETRICAL AND DESIGN PARAMETERS by Ehsan Zaman A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) October 2017 © Ehsan Zaman, 2017 ii Abstract This thesis describes a numerical study of the effect of swirl on the flow hydrodynamics in hydrocyclones and rotating pipes. Kinetic (force balance) and kinematic (velocity field) analysis were performed to identify the relative importance of turbulence, convective, and pressure forces in these flows. We showed that convective and pressure forces in hydrocyclones are three orders of magnitude larger than the size of turbulent forces, and seven orders of magnitude greater than the viscous forces. We derived a new dimensionless parameter based on a ratio of Euler numbers that gave a unique relationship for pressure along the axis of a hydrocyclone for all the cone angles tested. This finding enables rigorous scaling of hydrocyclones of differing cone angles. A complementary numerical study of flow in rotating pipes was carried out to elucidate the relative importance of convection and turbulence. We identified a dimensionless rotation parameter that delineates the condition at which decreasing turbulence force equals increasing convective force as rotational speed increases. This dimensionless number establishes a criterion for knowing which forces are dominant, and thereby a rational basis for choosing CFD (Computational Fluid Dynamics) models that are both cost-effective and accurate. Lastly, we conducted a sensitivity study to determine the effects of varying values of constants in the quadratic pressure-strain (QPS) Reynolds Stress turbulence model (RSM) for CFD modelling of swirl in hydrocyclones. The results identified the changes necessary to attain accurate predictions, such as overcoming the under-prediction problems of RSM. In addition, these findings laid the groundwork for a swirl-specific turbulence model for hydrocyclones. iii Lay Summary This thesis was written based on a numerical study of turbulent, swirling flows encountered in many industrial applications including separation devices such as hydrocyclones. The results provide a foundation for improving the computational modelling of fluid flows in swirling separation and classification equipment including hydrocyclones. A new dimensionless group was identified which facilitates the scaling of hydrocyclones, therefore creating a cost-effective design tool. Moreover, a novel criterion was proposed for establishing dominant forces in swirling, turbulent flows, thereby guidelines were provided for choice of turbulence model based on the proposed criterion. This provides engineers and researchers with a threshold to employ accurate models while reducing the computational cost in the process. In addition, by performing a sensitivity study we laid a foundation for fine-tuning the Reynolds Stress Model (RSM) into a hydrocyclone-specific turbulence model. Recommendations on favourable numerical values of the RSM model were made. iv Preface This thesis is an unpublished original work by Ehsan Zaman hereafter referred to as the “author”. This research project was proposed by Dr. James Olson and Dr. Mark Martinez. The research presented in this work was performed by the author. The author ran simulations and post-processed the computational fluid dynamics (CFD) data. Dr. Ali Vakil, a research colleague, helped analyze the numerical results. The experimental results are those of Dr. MacKenzie, as a research colleague at the Pulp and Paper Centre. v Table of Contents Abstract .......................................................................................................................................... ii Lay Summary ............................................................................................................................... iii Preface ........................................................................................................................................... iv Table of Contents ...........................................................................................................................v List of Tables ................................................................................................................................ xi List of Figures .............................................................................................................................. xii Nomenclatures ............................................................................................................................ xix List of Abbreviations ................................................................................................................ xxii Acknowledgements .................................................................................................................. xxiii Dedication ................................................................................................................................. xxiv Chapter 1: Introduction ................................................................................................................1 1.1 Historical Background .................................................................................................... 1 1.2 Applications .................................................................................................................... 3 1.2.1 Separation Based on Size ............................................................................................ 3 1.2.2 Separation Based on Density ...................................................................................... 3 1.3 Hydrocyclone Flow Field and Principles ........................................................................ 4 1.3.1 Tangential Velocity ..................................................................................................... 6 1.3.2 Flow Reversal ............................................................................................................. 9 1.3.3 Radial Velocity ......................................................................................................... 12 1.3.4 Short Circuit .............................................................................................................. 16 1.3.5 Mantle ....................................................................................................................... 17 vi 1.3.6 Particle Separation .................................................................................................... 19 1.3.7 Geometry Modifications and Effect of Design Parameters ...................................... 21 1.3.7.1 Effect of Vortex Finder Diameter ..................................................................... 25 1.3.7.2 Effect of Cone Length and Angle ..................................................................... 26 1.3.7.3 Effect of Inlet .................................................................................................... 26 1.4 Turbulence Modelling ................................................................................................... 27 1.4.1 Performance of Turbulence Models.......................................................................... 28 1.4.2 Effect of Swirl on Turbulence ................................................................................... 34 1.5 Summary of Literature and Objectives ......................................................................... 37 Chapter 2: Methodology..............................................................................................................39 2.1 Governing Equations .................................................................................................... 39 2.2 Turbulence Modeling .................................................................................................... 42 2.2.1 Standard k-ε Model ................................................................................................... 43 2.2.1.1 Reynolds Stress Tensor ..................................................................................... 44 2.2.1.2 Turbulent Transport and Pressure Diffusion..................................................... 45 2.2.2 Reynolds Stress Model ............................................................................................. 46 2.2.2.1 Turbulent Diffusion Transport .......................................................................... 47 2.2.2.2 Pressure-Strain Term ........................................................................................ 48 2.2.2.2.1 Linear Pressure-Strain ................................................................................. 48 2.2.2.2.2 Quadratic Pressure-Strain ............................................................................ 51 2.2.2.3 Dissipation Rate ................................................................................................ 52 2.2.3 Large Eddy Simulation (LES) .................................................................................. 52 2.2.3.1 Filtered Navier-Stokes Equations ..................................................................... 53 vii 2.2.3.2 Smagorinsky-Lilly SGS Model......................................................................... 57 2.2.3.3 Dynamic Model ................................................................................................ 58 2.3 Three-Dimensional Modeling ....................................................................................... 59 2.4 Axisymmetric Modeling ............................................................................................... 61 2.5 Ideal Vortex Laws ......................................................................................................... 62 2.5.1 Forced Vortex ........................................................................................................... 63 2.5.2 Free Vortex ............................................................................................................... 64 Chapter 3: Effect of Swirl on Turbulence .................................................................................65 3.1 Numerical Scheme and Boundary Conditions .............................................................. 66 3.2 Axisymmetric vs. Three-Dimensional Simulations ...................................................... 67 3.3 Grid Sensitivity ............................................................................................................. 68 3.4 Experimental Validation ............................................................................................... 69 3.5 Pressure Field ................................................................................................................ 70 3.6 Velocity Field................................................................................................................ 71 3.7 Force Balance Analysis................................................................................................. 73 3.8 An Integral Measure for Evaluating Convective forces vs Turbulence forces ............. 78 3.9 Guidelines on Choice of Turbulence Model for Swirling Flows .................................. 80 3.10 Summary ....................................................................................................................... 82 Chapter 4: Effects of Geometrical and Flow Parameters ........................................................84 4.1 Numerical Scheme and Boundary Conditions .............................................................. 86 4.2 Experimental Validation ............................................................................................... 86 4.3 Effect of Cone Angle .................................................................................................... 91 4.3.1 On Tangential Velocity ............................................................................................. 91 viii 4.3.2 On Pressure ............................................................................................................... 93 4.3.3 On Viscous Dissipation............................................................................................. 97 4.3.4 On Vorticity .............................................................................................................. 98 4.3.5 On G-Force ............................................................................................................. 100 4.4 Force Balance Analysis............................................................................................... 102 4.4.1 Vortex Force ........................................................................................................... 104 4.4.2 Bernoulli Gradient .................................................................................................. 106 4.5 Summary ..................................................................................................................... 108 Chapter 5: Effects of Vortex Finder Diameter, Number of Inlets and Inlet Positioning ....110 5.1 Effects of Vortex Finder Diameter ............................................................................. 110 5.1.1 On Tangential Velocity ........................................................................................... 110 5.1.2 On Pressure ............................................................................................................. 111 5.1.3 On Vorticity ............................................................................................................ 112 5.1.4 On Centrifugal Acceleration ................................................................................... 118 5.1.5 On Bernoulli Gradient............................................................................................. 120 5.2 Effects of the Number of Inlets and Inlet Positioning ................................................ 122 5.2.1 On Pressure ............................................................................................................. 124 5.2.2 On Tangential Velocity ........................................................................................... 125 5.2.3 On Centrifugal Acceleration ................................................................................... 126 5.3 Summary ..................................................................................................................... 127 Chapter 6: Sensitivity Study on the Reynolds Stress Model for the Hydrocyclone Flow Field .............................................................................................................................................128 6.1 Axisymmetric Assumption ......................................................................................... 128 ix 6.2 Computational Domain and Numerical Scheme......................................................... 128 6.3 Grid Sensitivity ........................................................................................................... 131 6.4 Experimental Validation ............................................................................................. 131 6.5 Effect of RSM Coefficients ........................................................................................ 134 6.5.1 Doubled Cases ........................................................................................................ 135 6.5.1.1 Insensitive ....................................................................................................... 137 6.5.1.2 Moderately Sensitive ...................................................................................... 137 6.5.1.3 Highly Sensitive .............................................................................................. 137 6.5.2 Halved Cases ........................................................................................................... 138 6.5.2.1 Insensitive ....................................................................................................... 140 6.5.2.2 Moderately Sensitive ...................................................................................... 140 6.5.2.3 Highly Sensitive .............................................................................................. 141 6.6 Hydrodynamics of the Flow Field .............................................................................. 141 6.7 Force Balance Analysis for LPS and QPS .................................................................. 143 6.8 Summary ..................................................................................................................... 145 Chapter 7: Summary and Conclusions ....................................................................................147 7.1 Objective 1 .................................................................................................................. 147 7.2 Objective 2 .................................................................................................................. 147 7.3 Objective 3 .................................................................................................................. 147 Chapter 8: Implications and Limitations of the Research and Recommendations for Future Work............................................................................................................................................149 8.1 Objective 1 .................................................................................................................. 149 8.2 Objective 2 .................................................................................................................. 149 x 8.3 Objective 3 .................................................................................................................. 149 References ...................................................................................................................................151 xi List of Tables Table 4.1: Hydrocyclones dimensions in mm............................................................................... 84 Table 4.2: Volume-averaged values of the dissipation function (ΦUT for r=3 mm) ..................... 98 Table 6.1: Hydrocyclone dimensions in mm. Labels have been illustrated in Figure 6.1. ......... 129 Table 6.2: RSM coefficients used in doubled cases; those highlighted in green are double the default values in case 0. .............................................................................................................. 135 Table 6.3: RSM coefficients used in halved cases; those in green are half the default values in case 0. .......................................................................................................................................... 138 xii List of Figures Figure 1.1: Schematic of hydrocyclone flow field structure (Habibian et al. 2008). ..................... 5 Figure 1.2: Schematic of hydrocyclone flow field and suspended particles (Cullivan et al. 2004).......................................................................................................................................................... 6 Figure 1.3: Schematic of short circuit and eddy flow (mantle), adopted from Bradley (1965). ... 16 Figure 1.4: Simplified schematic of the axial and radial flow pattern (Svarovsky, 1984). .......... 18 Figure 1.5: Hydrocyclone axial and radial flow pattern, visualized with tracer injection (Ohtake et al. 1987) .................................................................................................................................... 19 Figure 1.6: Measured tangential velocity distribution of seeding aluminum particles (Kelsall, 1952). .............................................................................................................................................. 8 Figure 1.7: LDA-measured tangential velocity distribution across the hydrocyclone (Dabir & Petty, 1986) ..................................................................................................................................... 9 Figure 1.8: Axial velocity distribution in a hydrocyclone (Kelsall, 1952); locus of zero vertical velocity in red. .............................................................................................................................. 10 Figure 1.9: Axial velocity distribution across a hydrocyclone (Dabir & Petty, 1986). ................ 11 Figure 1.10: Axial velocity distribution across a hydrocyclone measured using LDV (Monredon et al., 1992). .................................................................................................................................. 12 Figure 1.11: Radial velocity distribution across a hydrocyclone, adapted from Kelsall (1952). . 13 Figure 1.12: Radial velocity distribution inside a hydrocyclone; (left) with air core and (right) without air core adapted from Quian et al. (1989). ....................................................................... 14 Figure 1.13: (left) Radial velocity profiles, normalized by the average axial (through-flow) velocity at each plane and (right) schematic of the hydrocyclone and position of the investigated planes adopted from Fisher and Flack (2002)............................................................................... 15 xiii Figure 1.14: Regions of similar particles size distributions in a hydrocyclone adopted from Renner and Cohen (1978) (http://www.tandfonline.com). ........................................................... 20 Figure 1.15: Schematics of hydrocyclones used by Chiné and Concha (2000) with included cone angles of (a) β=14° and (b) β=154°. ............................................................................................. 22 Figure 3.1: Dimensionless axial velocity at x/D=20 at Re=7000 and N=0.29; normalized by the bulk axial velocity. ........................................................................................................................ 67 Figure 3.2: Dimensionless swirl velocity at x/D=20 at Re=7000 and N=0.29; normalized by the pipe wall circumferential velocity. ............................................................................................... 68 Figure 3.3: Axial velocity profiles for three different grids at x / D =20 at Re = 7000 and N = 0.69 Figure 3.4: Comparison between the normalized axial velocity profiles from the present work (PW), Hirai et al.’s RSM results, and Kikuyama et al.’s experimental data at Re = 10000 for two rotation numbers N = 0.5 and N = 1. ............................................................................................. 70 Figure 3.5: Euler number distribution along the dimensionless pipe length, normalized by diameter, at Re=7000 and a range of rotation numbers. ............................................................... 71 Figure 3.6: Axial velocity profiles, normalized by mean axial velocity, at different axial locations for Re=7 000 and N=1.43. ............................................................................................................. 72 Figure 3.7: Swirl velocity profiles, normalized by the pipe wall rotational velocity, at different axial locations for Re=7 000 and N = 1.43. .................................................................................. 72 Figure 3.8: Reynolds shear stress at Re = 7 000 and x = 19D for different rotation numbers. ..... 73 Figure 3.9: Accelerations, different forces in equation (2.73) normalized by density, for the stationary wall case (N=0) at different Reynolds numbers. .......................................................... 74 Figure 3.10: Non-dimensionalized accelerations, different forces in equation (2.73) normalized by the rotating wall centrifugal force, for low to intermediate rotation numbers at Re = 7 000. . 76 xiv Figure 3.11: Non-dimensionalized accelerations, different forces in equation (2.73) normalized by the rotating wall centrifugal force, for intermediate to high rotation at Re = 7000. ................ 77 Figure 3.12: All forces profiles, normalized by their respective pipe wall centrifugal force, for different Reynolds numbers at N=0.29. ........................................................................................ 77 Figure 3.13: Non-dimensionalized integral measures of the convective and turbulence forces, normalized by the integral measure for the pressure force, as a function of wall rotational speed for different Reynolds numbers. ................................................................................................... 79 Figure 3.14: Non-dimensionalized integral measures of the convective and turbulence forces, normalized by the integral measure of the pressure force, as a function of rotation number (N=Rω/Um) for different Reynolds numbers. The transition point between the forces ratios occurs at the same rotation number of Ntr ≈ 0.45. .................................................................................... 80 Figure 3.15: Comparison between the axial velocity predictions of k-ε model and RSM at Re=7000 for a range of rotation numbers, 0 ≤ N ≤2.86. ............................................................... 81 Figure 4.1: (a) The hydrocyclone general geometry; (b) overlay of the three cone angle geometries, cases A (red), B (green), and C (blue). ...................................................................... 85 Figure 4.2: Tetrahedral grid of case A. Similar grid was used for cases B and C. ....................... 85 Figure 4.3: Clockwise from top left; tangential velocity profiles at z1=71.9, z2=91, z3=110.1, and z4=129.1 mm below the inlet vs. the radial position normalized by the local wall radius. .......... 87 Figure 4.4: Clockwise from top left; tangential velocity profiles at z5=148.2, z6=167.3, z7=186.4, and z8=205.5 mm below the inlet vs. the radial position normalized by the local wall radius. .... 88 Figure 4.5: Clockwise from top left; tangential velocity profiles at z=71.9, z=91, z=110.1, and z=129.1 mm below the inlet vs. the radial position normalized by the local wall radius. ............ 89 Figure 4.6: Axial velocity; theoretical solutions vs. numerical solutions. .................................... 90 xv Figure 4.7: Vorticity magnitude contours at (a) mid-plane x=0, and (b) plane z1 for case A (1/s)........................................................................................................................................................ 91 Figure 4.8: Tangential velocity profiles at (a) z=71.9, (b) z=91, (c) z=110.1, and (d) z=129.1 mm below the inlet vs. the radial position normalized by the local wall radius. ................................. 92 Figure 4.9: Relative difference in the tangential velocity profiles of the cases C (θ = 15) and A (θ = 6.76) at z1 and z3. ....................................................................................................................... 93 Figure 4.10: Profiles of the pressure drop at z3. Radial position is normalized by the local wall radius. ............................................................................................................................................ 94 Figure 4.11: Pressure drop at z3 normalized by the pressure difference between underflow and overflow. ....................................................................................................................................... 95 Figure 4.12: The distribution of the ratio of the local Euler number to the global Euler number at z3. ................................................................................................................................................... 96 Figure 4.13: Pressure difference ratio (PDR) for the three cone angle cases. .............................. 96 Figure 4.14: Vorticity magnitude profiles at z3 vs. dimensionless radial position. ...................... 99 Figure 4.15: Radial distribution of centrifugal field for the three different cone angle cases at z3...................................................................................................................................................... 102 Figure 4.16: (right) Bernoulli gradient, vortex force, viscous force, and turbulence; (left) viscous and turbulence forces, at z3 for case A. The forces represent the terms in equation (4.9) multiplied by density, ρ=998.2 kg/m3. ........................................................................................ 103 Figure 4.17: A comparison between the two components of the vortex force at z3 for case A. . 104 Figure 4.18: Comparison between the axial vorticity terms of equation (4.11) at z3 for case A...................................................................................................................................................... 105 xvi Figure 4.19: The radial distribution of the kinetic energy gradient terms, LHS of equation (4.15), at z3 for case A. ........................................................................................................................... 107 Figure 4.20: The radial distribution of Bernoulli gradient at z3 for all cone angle cases. .......... 108 Figure 5.1: Tangential velocity profiles for the three vortex finder diameter cases at axial distances of (top left) z1=71.9 mm, (top right) z2=91.0 mm, (bottom left) z3=110.1 mm, and (bottom right) z4=129.1 mm below inlet. Radial position is normalized by the respective local wall radius. .................................................................................................................................. 111 Figure 5.2: Pressure Difference Ratio vs. the ratio of overflow diameter to underflow diameter for cases A, B, and C. ................................................................................................................. 112 Figure 5.3: Z-vorticity magnitude at z1=71.9 mm below the inlet, for vortex diameter cases. .. 114 Figure 5.4: Z-vorticity magnitude at z3=110.1 mm below the inlet (top of the cone) for vortex diameter cases. ............................................................................................................................ 114 Figure 5.5: Z-vorticity magnitude at z10=377.3 mm below the inlet, for vortex diameter cases. 115 Figure 5.6: Z-vorticity magnitude at various heights for case B. ............................................... 116 Figure 5.7: Z-vorticity magnitude at various heights for case A. ............................................... 116 Figure 5.8: Z-vorticity magnitude at various heights for case C. ............................................... 117 Figure 5.9: Centrifugal acceleration vs. normalized radius at z1 for the vortex finder diameter cases. ........................................................................................................................................... 118 Figure 5.10: Centrifugal acceleration vs. normalized radius at z4 for the vortex finder diameter cases. ........................................................................................................................................... 119 Figure 5.11: Centrifugal acceleration vs. normalized radius at z8 for the vortex finder diameter cases. ........................................................................................................................................... 119 xvii Figure 5.12: Centrifugal acceleration vs. normalized radius at z9 for the vortex finder diameter cases. ........................................................................................................................................... 120 Figure 5.13: Bernoulli gradient vs. normalized radial position at z1 for cases A, B, and C. ...... 121 Figure 5.14: Bernoulli gradient vs. normalized radial position at z8 for cases A, B, and C. ...... 122 Figure 5.15: Bernoulli gradient vs. normalized radial position at z9 for cases A, B, and C. ...... 122 Figure 5.16: Geometries of hydrocyclone (a) case D with s=0 and (b) case E with s=d. ........... 123 Figure 5.17: Pressure drop profiles at heights z1, z8, z9, z10; radial position normalized by the local wall radius. ......................................................................................................................... 124 Figure 5.18: Normalized pressure drop at heights z1, z8, z9, z10; radial position normalized by the local wall radius. ......................................................................................................................... 125 Figure 5.19: Tangential velocity profiles at heights, clockwise from top left, z1, z8, z9, z10; radial position normalized by the local wall radius. ............................................................................. 126 Figure 5.20: Centrifugal acceleration profiles at height z1 vs. normalized radial position. ........ 127 Figure 6.1: The computational domain and mesh used in the simulations, mirrored with respect to the axis. The high concentration of mesh elements makes it difficult to clearly see the grid structure; thus, the zoomed-in cells are displayed for the sake of clarity. .................................. 129 Figure 6.2: Swirl velocity profile for different grid sizes at z = 364.5 mm. ............................... 131 Figure 6.3: Swirl velocity profiles at (a) z=97.1 mm, (b) z= 115.4 mm, (c) z=135.5 mm, and (d) z=154.5 mm, away from the cyclone roof; radial position is normalized by the local radius. ... 132 Figure 6.4: Axial velocity profiles at (a) z=97.1 mm, (b) z= 115.4 mm, (c) z=135.5 mm, and (d) z=154.5 mm, away from the hydrocyclone roof; radial position is normalized by the local radius...................................................................................................................................................... 133 xviii Figure 6.5: Swirl velocity profiles at (a) z=211.8 mm, (b) z= 250 mm, (c) z=370.3 mm, and (d) z=364.5 mm, away from the cyclone roof; radial position is normalized by the local radius. ... 134 Figure 6.6: Swirl velocity profiles at z=85.4 mm below the roof (in the cylinder) for the QPS cases with increased constants vs. the default case..................................................................... 136 Figure 6.7: Swirl velocity profiles at z=225.40 mm below the roof (in the cone) for the QPS cases with increased constants vs. the default case..................................................................... 136 Figure 6.8: Swirl velocity profiles at z=85.4 mm below the roof for the QPS cases with reduced constants vs. the default case. ..................................................................................................... 139 Figure 6.9: Swirl velocity profiles at z=225.4 mm below the roof for QPS cases with reduced constants vs. the default case. ..................................................................................................... 140 Figure 6.10: Profiles of all forces in equation (2.74) (left), and profiles of viscous and turbulence forces (right) at z=102 mm from the cyclone roof. These results correspond to the QPS model...................................................................................................................................................... 142 Figure 6.11: Pressure and convective forces profiles for the QPS and LPS models at z=102 mm away from the roof - the border of the conical and cylindrical bodies. ...................................... 144 Figure 6.12: The ratio of the centrifugal force to the convective forces in the QPS model, at z=102 mm from the roof. ............................................................................................................ 144 xix Nomenclatures Roman Nomenclature D Diameter 𝑉𝑟 Radial velocity 𝑉𝜃 Tangential velocity 𝑉𝑧 Axial velocity r Radial position k Kinetic energy N Dimensionless rotation number R Wall radius 𝑈𝑚 Mean axial velocity x Axial distance L Pipe length ui Velocity component g Gravitational acceleration p Pressure 𝑝′ Fluctuating pressure P Mean pressure 𝑠𝑖𝑗 Strain-rate tensor 𝑠𝑖𝑗′ Fluctuating stress tensor 𝑈𝑖 Mean velocity component Τ𝑖𝑗 Mean stress tensor xx 𝑢𝑖′ Fluctuating velocity component 𝑏𝑖𝑗 Reynolds anisotropy tensor Eu Euler number Re Reynolds number a Acceleration di Inlet diameter do Overflow diameter du Underflow diameter Lv Vortex finder length Lcyl Cylindrical section length Lcon Conical section length t Time 𝑈𝑟 Radial velocity 𝑈𝜃 Tangential velocity 𝑈𝑧 Axial velocity B Bernoulli function Greek Nomenclature 𝜔 Angular velocity ε Dissipation rate ρ Density 𝜏𝑖𝑗 Stress tensor 𝜏𝑖𝑗′ Instantaneous stress tensor xxi 𝜏𝑖𝑗𝑡 Reynolds stress tensor 𝜈 Kinematic viscosity μ Dynamic viscosity 𝜇𝑡 Turbulent eddy viscosity Φ𝑖𝑗 Pressure-strain tensor Ω𝑖𝑗 Mean rotation-rate tensor δ Kronecker delta Δ Filter cutoff width θ Cone angle 𝜔𝑧 Vorticity in the z-direction ?⃗? Vorticity vector Ω𝑧 Mean z-vorticity xxii List of Abbreviations LZVV Locus of zero vertical velocity LDA Laser Doppler anemometry LDV Laser Doppler velocimetry CFD Computational fluid dynamics RSM Reynolds stress model LES Large eddy simulation RNG ReNormalization Group LPS Linear pressure-strain QPS Quadratic pressure-strain DNS Direct numerical simulations RANS Reynolds-Averaged Navier-Stokes LHS Left-hand side RHS Right-hand side PS Pressure-strain PDE Partial differential equation 2D Two-dimensional 3D Three-dimensional TVRD Tangential velocity relative difference PDR Pressure difference ratio xxiii Acknowledgements I would like to express my sincere gratitude to my supervisor Dr. James Olson and supervisory committee member Dr. Mark Martinez for all their support and patience during the course of this project. I am immensely grateful to you both. I would also like to sincerely thank Dr. Richard Kerekes for his assistance and invaluable feedback during this project. I would also like to thank the National Sciences and Engineering Research Council of Canada (NSERC) and FPInnovation for their financial support. The staff and students at the Pulp and Paper Centre (PPC) have been supportive throughout this process. Special thanks to George Soong for his help and devotion to myself and the students at PPC. I would also like to thank my friend and research colleague Dr. Ali Vakil for many helpful discussions and words of encouragement throughout the course of this project. Many thanks to my dear friend and colleague Hamed Ghasvari for his invaluable assistance and support over the years. I would like to extend my appreciation to other colleagues and friends at PPC, including but not limited to, Dr. Frank Saville, Dr. Abbas Nikbakht, Mohammad Shanbghazani, Justin Roberts and Hatef Rahmani. I owe a great deal of appreciation to my family for their unconditional love and support throughout my years of study. Special thanks to my dad, Malek, and my sisters Neda, Nasim, and Nahal. Most importantly, my mom has always been my inspiration. Thank you for all you did for me and I hope I have made you proud. I would like to pay special tribute to my loving partner, Vanessa, for her support, patience, and understanding throughout this project. I would also like to thank Vanessa for proofreading this thesis and putting up with me during periods of frustration. Lastly, I would like to thank my cute, adorable dogs, Alita and Arty, for being by my side while putting this thesis together. xxiv Dedication In memory of my mother. 1 Chapter 1: Introduction A comprehensive review of the literature on the hydrocyclone flow field is presented. Particular attention is paid to turbulence modelling, and the effect of the hydrocyclone design variables on its flow characteristics. Additionally, a general summary of studies on the effect of rotation (swirl) on turbulence is given. A thorough understanding of the interaction between swirl and turbulence will greatly help gain an in-depth understanding of turbulent, swirling flows, including the hydrocyclone flow field. 1.1 Historical Background Hydrocyclones are separation devices which have found wide industrial applications over the past several decades. The term hydrocyclone is an abbreviation for ‘hydraulic cyclone’, which is a reference to water being the fluid medium. Gas cyclones, gas being the fluid medium, had been in use for several decades prior to the hydrocyclone. The basic principle behind cyclone separation is centrifugal sedimentation; particles subjected to the centrifugal field are separated or classified in the fluid domain. The same principle for separation is used in the centrifuges; however, there is a major difference when it comes to cyclones in that the required vertical motion is generated by the fluid itself as cyclones have no moving parts. The hydrocyclone was patented over a century ago (Bretney, 1891). A U.S. patent granted to MacNaughton (1906) was the first time that the hydrocyclone was applied to the pulp and paper industry. That was followed by a British patent (Berges, 1935) for the use of hydrocyclones for cleaning the pulp stock. Shortly afterwards, Freeman (1936) was granted another U.S. patent for the work developed at Consolidated Papers in Quebec. Several patents were obtained by the Dutch State Mines in the Netherlands for the use of hydrocyclones in the mining and mineral processing industries including the dewatering of fine sand suspension –which would act as the heavy medium for coal-shale separation (Bradley 1965). Immediate success of the hydrocyclone use in the mineral industries led to a broader range of applications for the hydrocyclone. Industries including the pulp and paper, textile, chemical and petroleum, adopted the hydrocyclone for separation and classifying purposes among others. The first English text on the hydrocyclone subject was by Bradley (1965). A few decades later, a “revival of interest” in hydrocyclones motivated Svarovsky 2 (1984) to write a full text on the topic. The revived interest was attributed to factors such as interdisciplinary communication between mineral processing and other fields including chemical engineering and process engineering. Additionally, other factors included the need of the petroleum industry for separating sand, gas, or water from oil using simple and robust tools, and the increase of the overall knowledge of the hydrocyclone flow pattern. The advantages of using hydrocyclones include (Svarovsky 1984), I. They are simple and stationary with no moving parts. II. They are very cost-effective both in terms of capital investment and installation as well as operation and maintenance. III. They are relatively small in size, thus less space-restrictive compared to other separation devices. IV. The high shear forces existing in the flow breaks agglomerates, thus providing an advantage to the hydrocyclone in classification of solids, and in treatment of Bingham plastic slurries. V. Their relatively low residence times is an advantage over sedimentation classifiers in terms of the speed of control. VI. They are exceptionally versatile in application, which makes them attractive to a broad range of industries and with a wide variety of uses. They are used for classifying solids, washing solids, clarifying liquids, and separating immiscible liquids among others. Subsection 1.2 provides additional details about these applications. The main disadvantages of using hydrocyclones can be summarized as follows (Svarovsky 1984), I. Their separation performance may be somewhat poor in terms of sharpness and range of operating cut size. Multistage configurations may help overcome some of these shortcomings, but that will add to the costs. II. They are rather inflexible devices due to the influence of the flow rate and feed concentration on their separation performance. III. They are abrasion prone; however, measures can be taken to minimize the impact. 3 IV. Although mostly an advantage, the presence of high shear forces can become a disadvantage; flocculation cannot be used unlike gravity thickeners, which exploit flocculation to enhance the separation performance. 1.2 Applications Hydrocyclones are employed in various fields of technology and a wide range of industries with numerous applications including clarification, thickening, and solid-solid separation. The results of this work will be applicable chiefly to solid-solid separation with low particles concentration. Hence, the following subsection elaborates on this particular application. See Svarovsky (1984) for more details on hydrocyclone applications. 1.2.1 Separation Based on Size The aim of solid-solid separation according to size is to fractionate solid particles into coarse and fine streams. Alternatively and more commonly, this form of separation is desirable to either remove coarse particles, such as in refining operations, or to remove fine particles, for instance in washing operations. This separation process can be done in multiple stages in an attempt to compensate for the poor sharpness in the hydrocyclone fractionated streams. Additional examples of this category include, de-gritting of drilling mud and de-gritting of underground waters classifying sand according to size classifying iron ore into sinter stream (coarse) and pellet stream (fine) 1.2.2 Separation Based on Density Solid-solid separation according to density separates solids based on their density. This method is common in mining and mineral processing industries for sorting minerals, coal and rocks by either the dense medium cyclone or water-only cyclones. It is worth mentioning that the solid-solid separation, in certain applications, can be achieved based on both the particle shape and density. The cleaning of pulp slurry for papermaking is the 4 most widely used example of solid-solid separation according to shape and density. Other examples include, fibre fractionation in pulp suspensions removal of wax, sand, clay, bark, shives, and other contaminants from pulp slurry Immiscible liquids can be separated successfully using cyclones, for example the separation of oil from water and dewatering light oils. Another application of cyclones is degassing of crude oil, and in general gas bubbles separation from liquids. A thorough understanding of the flow field inside the hydrocyclone is crucial for fully comprehending its separation/classification mechanism, and for optimal design of the hydrocyclone intended for a specific application. A chronological literature summary of the hydrocyclone flow field is given in the following subsection. Certain aspects of its hydrodynamics which are not clearly understood yet, require additional investigation to enable engineers to further advance and optimize the design for desired outcome. These aspects are discussed below. 1.3 Hydrocyclone Flow Field and Principles The schematic of the design and flow field pattern of the conventional, cono-cylindrical (or frusto-conical) hydrocyclone has been depicted in Figure 1.1 and Figure 1.2). This type of hydrocyclone, the most widely used in industry and research, is composed of a short cylindrical body followed by a long cone, with the outlets at opposite ends, usually referred to as overflow (top exit) and underflow (bottom exit). Some of the fluid flow exits the hydrocyclone via the overflow through the vortex finder. The vortex finder is the opening at the end of a pipe protruding inside the hydrocyclone through the top wall (sometimes called roof) of the cylinder. The rest of the fluid flow leaves the hydrocyclone through the underflow orifice, which is the opening at the apex of the cone. Fluid enters the cylinder through a tangential inlet, typically either a circular or rectangular pipe. As a result of the tangential entry, a strong swirling motion is generated within the hydrocyclone and the pressure head of the feed stream is converted into rotational motion. Consequently, as shown in Figure 1.1, a downward helical flow pattern (the outer helical or outer vortex) formed in the wall proximity, exits through the underflow, while an upward helical flow 5 pattern (the inner helical or vortex flow) formed near the hydrocyclone central axis, leaves through the vortex finder tube and subsequently the overflow opening. The outer helical carries larger and heavier particles. Figure 1.1: Schematic of hydrocyclone flow field structure (Habibian et al. 2008). 6 Figure 1.2: Schematic of hydrocyclone flow field and suspended particles (Cullivan et al. 2004). 1.3.1 Tangential Velocity Tangential velocity is the most dominant velocity component in the hydrocyclone flow field. It has been investigated most widely, compared to the axial and radial components, mainly due to its strong influence on particle separation, and partly due to the more readily available techniques to measure and analyze the tangential velocity field. Virtually, all the initial studies on the hydrocyclone flow field were performed using experimental measurement methods. One of the first and most established studies in the literature was conducted by Kelsall (1952). He used fine aluminum particles to seed the water flow inside a 76 mm hydrocyclone, with the split ratio changing from 0 (the entire flow exiting through the underflow) to ∞ (the entire flow exiting through the overflow). The tangential velocity of the particles was directly measured using an optical microscope. Kelsall found that the tangential velocity increased with radial position, from the geometrical axis towards the wall, and reached its peak at a certain radius, before a sharp transition in its trend. The tangential velocity declined with further increases in radial position for radii greater than that corresponding to the peak. Kelsall’s measurements indicated the co-existence of a forced vortex and a free vortex. The inner vortex, within which the swirl velocity 7 was linearly proportional to the radius, closely resembled the solid-body motion, whereas the outer vortex exhibited characteristics of a free-vortex flow pattern. The tangential velocity pattern could be mathematically written as, 𝑉𝜃 = {𝑟𝜔, 𝑟 < 𝑟𝑝𝑒𝑎𝑘𝐶𝑟𝑛, 𝑟 > 𝑟𝑝𝑒𝑎𝑘 (1.1) where 𝑉𝜃 and r are tangential velocity and radial position, C is a constant, 𝜔 is angular velocity, n is an exponent, and subscript peak indicates the radial position of the tangential velocity peak. 𝑟𝑝𝑒𝑎𝑘 was found to be smaller, at all heights, than the radius of the vortex finder tube. The numerical value of the exponent n suggested by Kelsall, was between 0.75 and 0.84 for the split ratio between 0 and ∞, respectively. Kelsall observed that the tangential velocity was almost independent of axial height; envelopes of constant tangential velocity form co-axial cylinders around the cyclone geometrical axis, as illustrated in Figure 1.3. Yoshioka and Hotta (1955) reported a value of 0.8 for n for a dilute slurry flow inside a hydrocyclone. Ohasi and Maeda (1958) used a hydrocyclone of the same size as Kelsall’s (but of different design) and reported an average value of 0.74 for n. Knowles et al. (1973) used high-speed cine photography to measure the velocity components of tracer anisole droplets of 50-300 μm in size and a density of 993 kg/m3. While the qualitative trend of the tangential velocity was similar to that of Kelsall’s, the numerical value of the exponent differed significantly; n varied between 0.2 and 0.4 compared to 0.75 and 0.84 of Kelsall’s. 8 Figure 1.3: Measured tangential velocity distribution of seeding aluminum particles (Kelsall, 1952). Dabir and Petty (1986) employed Laser Doppler Anemometry (LDA) to measure the velocity distribution in a hydrocyclone with similar design and operating conditions to Knowles et al. (1973), and a split ratio of 4, i.e. 80% of the flow exited through the overflow. Their results showed that varying the axial position led to little changes in the tangential velocity profiles. The tangential velocity profile remained almost unchanged across the hydrocyclone irrespective of axial position, as shown in Figure 1.4. Hence, the tangential velocity profile was independent of axial position, which is consistent with Kelsall’s observations. The numerical value of exponent n was 0.62, which was substantially different than that of Knowles et al.(1973), despite similar conditions. 9 Figure 1.4: LDA-measured tangential velocity distribution across the hydrocyclone (Dabir & Petty, 1986) Other well established research works have yielded similar observations (Quian et al., 1989; Monredon et al., 1992; Chiné and Concha, 2000; Fisher and Flack, 2002; Concha, 2007). Most of these studies were conducted in pure or dilute suspensions. The presence of solid particles is expected to lower the maximum tangential velocity and leads to a smoother transition between the free vortex and forced vortex flows. This was shown to be the case for pulp suspensions at low fibre concentration (Bergstrom et al., 2007). 1.3.2 Flow Reversal The flow reverses its direction, from downward near the wall to upward near the centre, at a certain radial position inside the hydrocyclone. This flow reversal was shown by Kelsall (1952) and can be seen in Figure 1.5. The flow is downward in the outer vortex while it is upward in the inner vortex. Kelsall identified a locus of zero vertical velocity (LZVV) between these two counter-current vortices, which roughly speaking followed the hydrocyclone outline. The LZVV has been highlighted by dashed red lines in Figure 1.5. 10 Figure 1.5: Axial velocity distribution in a hydrocyclone (Kelsall, 1952); locus of zero vertical velocity in red. Knowles’ measurements supported Kelsall’s observation. However, Knowles et al. (1973) reported substantial fluctuations for the axial velocity inside the vortex finder tube; nevertheless the flow direction did not change. Dabir and Petty (1986) discovered that the axial velocity reversed direction more than once along the radial direction, as shown in Figure 1.6. 11 Figure 1.6: Axial velocity distribution across a hydrocyclone (Dabir & Petty, 1986). Ohtake et al. (1987) reported the existence of a downward flow just outside the air core. The LDV (Laser Doppler Velocimetry) measurements of Monredon et al. (1992) revealed that the axial velocity field had an axisymmetric behavior within the conical section, with asymmetries observed particularly in the top half of the cylindrical section, as illustrated in Figure 1.7. They attributed the existence of the asymmetry in the upper half of the cylinder to the geometrical asymmetry introduced by the single inlet. 12 Figure 1.7: Axial velocity distribution across a hydrocyclone measured using LDV (Monredon et al., 1992). 1.3.3 Radial Velocity The radial velocity is the smallest velocity component in the hydrocyclone flow field. Research literature on the radial velocity field is not as rich as its axial and tangential counterparts. Kelsall (1952) obtained radial velocity distribution at various heights by solving axisymmetric continuity equation using his tangential and axial velocity measurements. The calculated radial velocity profiles, as shown in Figure 1.8, indicated an inward distribution which diminished with 13 decreasing radius. The findings of Knowles et al. (1973) was significantly different than Kelsall’s. They employed geometric calculations to obtain the radial velocity distribution. They reported that the radial velocity was outward near the centre. Moreover, the magnitude of the radial velocities of Knowles et al. (1973) was lower than that of Kelsall (1952). Figure 1.8: Radial velocity distribution across a hydrocyclone, adapted from Kelsall (1952). Quian et al. (1989) used LDA to measure all three velocity components in two hydrocyclones; an ordinary hydrocyclone (with air core) and a water-sealed hydrocyclone (without air core). Their 14 measurements indicated that the radial velocity was inversely proportional to radius, as shown in Figure 1.9. This observation was in ‘complete’ disagreement with that of Kelsall whose radial velocity was directly proportional to radius. Figure 1.9: Radial velocity distribution inside a hydrocyclone; (left) with air core and (right) without air core adapted from Quian et al. (1989). Moreover, they found that the radial velocity of the water-sealed hydrocyclone was, on average, 0.5 m/s higher than that of the ordinary one. Also, Quian et al. (1989) expressed their measured radial velocity (𝑉𝑟) profile as, 15 𝑉𝑟𝑟𝑚 = −𝐶 (1.2) where m and C were positive constants. They stated that despite the dependency of exponent m and constant C on vertical position, the equation (1.2) held as the overall mathematical form for the description of the radial velocity distribution. Note that the minus sign on the right-hand side of the equation indicates that the flow is inward, i.e. negative radial velocity. In a very interesting LDA study, Fisher and Flack (2002) measured all three velocity components for a variety of inlet flow rates and reject rates in a forward flow hydrocyclone. They reported that the magnitudes and shapes of the radial velocity profiles showed no significant dependency on axial position, except near the air core, as illustrated in Figure 1.10. Their measurements indicated that the radial velocity was nearly zero, at all axial heights, but near the air core. Figure 1.10: (left) Radial velocity profiles, normalized by the average axial (through-flow) velocity at each plane and (right) schematic of the hydrocyclone and position of the investigated planes adopted from Fisher and Flack (2002). The mean of the absolute values for radial velocities was claimed to be only 2% of the mean of the tangential velocities. They concluded that the radial velocity was the least dominant velocity component in terms of magnitude, however, it was the most critical for the separation process. 16 The flow field inside the hydrocyclone is categorically ‘quasi-steady’ or ‘quasi-periodic’ rather than steady. This is due to the periodical fluctuations constantly generated by the precession of the vortex core (PVC). The PVC is suspected to influence the particle separation efficiency. However, the extent of its influence is not clear (Bergström, 2006). 1.3.4 Short Circuit Short circuit flow is a secondary flow pattern at the top of the cylindrical body which provides a short and direct flow path from the inlet to the vortex finder tube, across the cyclone roof (top wall) and down the outer wall of the vortex finder. The presence of the roof, and the vortex finder outer wall, leads to a local slow-down of the swirl of the feed stream. This slowdown creates a low-resistance flow area, adjacent to the roof and outer wall of the vortex finder, for flow from the outer regions with far higher pressure. The coupled effects of the locally reduced swirl and higher pressure in the outer regions, forces some portion of the feed stream to pass directly to the vortex finder, across the cyclone roof and down the vortex finder outer wall. A simplified illustration of short circuit flow has been provided in Figure 1.11 -Figure 1.13). Short circuit is undesirable, as particles join the overflow through this direct path without experiencing the hydrocyclone centrifugal field, and without being classified or separated. Figure 1.11: Schematic of short circuit and eddy flow (mantle), adopted from Bradley (1965). 17 1.3.5 Mantle The entirety of the upward fluid flow in the vicinity of the vortex finder opening does not manage to pass through the tube. Subsequently, the portion of the upward flow not entering the vortex finder tube continues its upward motion just outside the outer wall of the vortex finder. At some point the flow direction of motion changes from upward to downward; therefore, this flow pattern forms one or more recirculating eddies, as depicted in Figure 1.11. These recirculation eddies (zones) are referred to as the mantle and are marked in Figure 1.12 andFigure 1.13. There is little-to-zero radial flow through the mantle. Hence, the mantle can be regarded as a buffer zone between the downward flow of the outer helical and the upward flow in the inner helical. Using the dye injection technique, Bradley and Pulling (1959) reported that increasing the cone angle yielded to stronger recirculation and the formation of additional eddies. They reported the position of the mantle, a cylindrical surface at a diameter of about 0.43D (the cyclone diameter) which always vanished at a height in the conical body corresponding to about 0.7D, were neither influenced by varying the cyclone diameter between 38 mm and 76 mm nor by the cyclone design within “normal” limits (Svarovsky 1984). Moreover, it was discovered that the radial flow appeared to exist along the entire length of the cyclone when the vortex finder tube radius exceeded the position of the mantle; thus, a stable, stationary mantle ceased to exist. 18 Figure 1.12: Simplified schematic of the axial and radial flow pattern (Svarovsky, 1984). 19 Figure 1.13: Hydrocyclone axial and radial flow pattern, visualized with tracer injection (Ohtake et al. 1987) 1.3.6 Particle Separation The classification or separation of particles (or fractionation) takes place as a direct consequence of the hydrocyclone flow field. Surprisingly, there has not been extensive research performed to directly relate the fluid flow field inside the hydrocyclone to the particle separation mechanism. Also, the distribution of particles within a hydrocyclone is far from known; often the focus of research has been put on the size distribution in the final products (accepts and rejects) rather than inside the hydrocyclone. Renner and Cohen (1978) conducted an experimental study to shed light on particle distribution inside a 150 mm diameter hydrocyclone. They developed a high-speed sampling probe; the experiment was designed so that the collected samples were unaffected by any disturbance caused 20 by the probe. This was achieved by careful design of the method and the speed of sampling. They discovered that the hydrocyclone can be divided into four zones of similar size distributions, as depicted in Figure 1.14. They identified zone A as a region which contained particles with the same size distribution as the inlet feed and particles in zone B had a very close size distribution to those in the underflow (rejects or coarse products). The particle size distribution in zone C was identical to that of the overflow (accepts or fine products), and toroidal-looking zone D was richer in intermediate size fractions than the rest of the hydrocyclone including the inlet feed. Put in terms of fractionation, zone A contained unfractionated particles, zone B contained fractionated coarse stream, zone C, within the vortex finder vicinity, contained fractionated fine stream, and finally zone D appeared to be the only region where fractionation occurred. Figure 1.14: Regions of similar particles size distributions in a hydrocyclone adopted from Renner and Cohen (1978) (http://www.tandfonline.com). 21 1.3.7 Geometry Modifications and Effect of Design Parameters A great deal of research has been performed over the years to understand the hydrocyclone flow field (Monredon et al., 1992; Chiné and Concha, 2000; Fisher and Flack, 2002; Nowakowski and Dyakowski, 2003; Cullivan et al., 2004; Petty and Parks, 2004; Roldan-Villasana et al., 1993). Chen et al. (2000) examined several existing, widely used theoretical and empirical hydrocyclone models. The study concluded that none of the models were able to accurately predict all the flow field characteristics. Additionally, no model was universal and their accuracy was limited to certain types of geometrical and operational features. Review studies of the hydrocyclone flow field modelling are available in the literature (Nowakowski et al., 2004; Narasimha et al., 2007). It has been established that the hydrocyclone geometry plays a critical role in its separation performance (Wang and Yu, 2006; Mousavian and Najafi, 2009; Vieira et al., 2011). Several studies have been conducted on the effect of geometry modifications so as to design hydrocyclones to meet certain separation demands in specific industries (Wang and Yu, 2008; Mokni et al., 2015; Delgadillo and Rajamani, 2007; Ghodrat et al., 2014a). Monredon et al. (1992) measured the velocity field inside five hydrocyclones by means of LDV. Their results indicated that the tangential velocity was higher in the design with larger cone angle and smaller vortex finder diameter. The difference between the tangential velocity profiles reduced at lower sections of the hydrocyclone; nevertheless, the tangential velocity of the larger cone angle and smaller overflow design remained the greatest across the hydrocyclone amongst the five investigated designs. The axial velocity field, however, showed a different trend. While the axial velocity profiles were similar at the top of the hydrocyclones, the axial velocity was highest in the design with smaller cone angle and apex diameter. Their results indicated that the design with smaller cone angle had the largest cut size for removal of limestone particles. Using LDV, Chiné and Concha (2000) measured and compared the axial and tangential velocity fields for water flow inside a 102 mm diameter hydrocyclone. The hydrocyclone body could be altered by switching to/from a conventional cono-cylindrical hydrocyclone (with a conical section) from/to a cylindrical (or flat bottom) hydrocyclone (with no conical section, also referred to as zero-degree cone angle), as illustrated in Figure 1.15. 22 Figure 1.15: Schematics of hydrocyclones used by Chiné and Concha (2000) with included cone angles of (a) β=14° and (b) β=154°. While the tangential velocity exhibited similar trends in the two types, the axial velocity patterns were drastically different. The axial velocity was a linear function of axial position in the cylindrical hydrocyclone, which was not the case for the conical hydrocyclone. Additionally, they found that turbulence intensity at the upper regions of the hydrocyclone was higher in the cylindrical hydrocyclone. However, that changed near the apex where turbulence intensity in the conical hydrocyclone exceeded that of the cylindrical vessel. Chu et al. (2000) experimentally investigated the effect of geometrical modifications and showed the importance of the conical body for the separation performance. They reported that a parabola cone type, as opposed to the standard straight cone, could increase the separation efficiency and sharpness as well as the cut size. Vieira et al. (2005) experimentally studied the separation performance of two different filtering hydrocyclone designs and highlighted the influence of the conical body on the hydrocyclone separation performance. Experimental studies on understanding 23 the effect of geometrical features are not economically viable because manufacturing new designs (or geometry modifications) is time-consuming and incurs high costs in experimentation. Alternatively, Computational Fluid Dynamics (CFD) offers a great tool, at a much lower cost, for optimizing the hydrocyclone geometry for specific applications. Wang and Yu (2006) used the Reynolds Stress Model (RSM) to study the effect of the hydrocyclone dimensions, including the length of the cylindrical and conical sections, on the hydrocyclone flow field. They claimed that a longer conical section could enhance the separation efficiency, whereas the length of the cylindrical section was reported to play an inconsequential role in the particle separation mechanism. The study pointed out that the hydrocyclone dimensions are crucial in the flow field pattern, as changes in geometrical sizes changed the flow field and separation performance. They concluded that CFD should be used a reliable tool to optimize the hydrocyclones dimensions for optimum performance in the application of interest. Delgadillo and Rajamani (2007) conducted numerical simulations on the effect of novel designs on the hydrocyclone separation efficiency. Their principal contribution was to vindicate CFD as the right tool for exploring new designs before experimental verification was performed. Additionally, they reported the predicted separation was sharper in a two-cone design, in which a 20-degree cone connected the cylindrical body to a second and longer conical section with an included angle of 6 degrees, than that of a standard hydrocyclone. They did not verify this result with experiments. In a separate CFD study, Wang and Yu (2008) focused on the effect of the shape and size of the vortex finder tube on the separation performance of a 75-mm diameter hydrocyclone. They observed opposite effects of the vortex finder tube length on the separation efficiency for fine and coarse particles. Due to the short circuit phenomenon, shortening the vortex finder enhanced the separation performance for coarser particles, and was detrimental to the separation efficiency for finer particles. It was also reported that in spite of weakening the eddy flow effects, thickening the vortex finder wall adversely affected the separation performance by strengthening the short circuit. Conversely, a thinner vortex finder wall led to improved separation performance; however, it caused a high pressure drop. 24 Hwang et al. (2008) conducted a CFD study on the effect of the hydrocyclone’s geometrical structure on the separation of fine particles. They reported that decreasing the vortex finder diameter and enlarging the underflow conduit could improve the separation efficiency. However, those changes would decrease the separation sharpness. Mousavian and Najafi (2009) numerically investigated the effect of cone angle, and both the underflow and overflow diameters on the hydrocyclone performance. The separation efficiency was reported to drop by decreasing the overflow diameter. They observed an enhanced separation performance for hydrocyclones with a longer conical section. Furthermore, they reported that the cut size decreased with an increase in cone angle or a decrease in the conical section, which was attributed to a stronger short circuit. The study concluded that the hydrocyclone flow field and separation performance is dependent on the dimension of different geometrical parameters. Using the RSM, Murthy and Bhaskar (2012) conducted a CFD parametric study on the flow field for a 76 mm diameter hydrocyclone with a 20-degree cone angle. They reported that widening the vortex finder opening led to a lower tangential velocity peak, a lower pressure differential, and a greater cut size. An increase in the underflow opening size was observed to have a similar effect, of a smaller scale, on the pressure and tangential velocity fields; however, it resulted in a lower cut size and a reduction in the separation efficiency. Also, increasing the throughput (flow rate) was reported to significantly increase the tangential velocity and pressure differential, which resulted in an improved separation efficiency and in a lower cut size. Considerable declines occurred in the tangential velocity and pressure differential when the fluid viscosity increased. Xu et al. (2013) performed a CFD study on the effect of the ratio of the overflow diameter to the underflow diameter on the air core formation and stability. They claimed that the diameter ratio was the decisive factor in the stability/formation of the air core, whereas operational conditions did not affect the steady state of the air core. They reported that at a ratio of 1.2 no air core was formed. Saidi et al. (2013) used LES (large eddy simulation) to compare the flow fields inside a 35 mm diameter de-oiling hydrocyclone with three different cone angles. They reported that larger cone 25 angles led to higher tangential velocity and pressure gradients, but decreased the separation efficiency due to reduced residence time. Hwang et al. (2013) employed RSM and investigated the effect of the number of inlets and their size on particle separation for the same flow rate. The inlet arrangements included single inlet, dual inlet (with half-size inlet and the same inlet velocity), and two tetra inlet cases. The first tetra-inlet case had identical inlet size as the original single-inlet case, meaning that the feed velocity had to be adjusted to retain the same flow rate. The size of the inlet in the second tetra-inlet case was reduced to the same feed velocity as the original case, which resulted in the same flow rate. Their study concluded that increasing the number of inlets with reducing the inlet size enhanced the separation efficiency. In a numerical study of the effect of the conical body length and shape on the hydrocyclone performance, Ghodrat et al. ( 2014a) used RSM and reported an improved separation performance for a hydrocyclone with a convex cone design, compared to the standard or a concave cone shape. The hydrocyclone separation performance can be enhanced through various means, such as increasing the throughput (flow rate), enlarging the conical section, reducing the hydrocyclone diameter. However, some changes could adversely affect other performances or requirements. For instance, increasing the flow rate would result in a higher pressure drop, hence a higher loss. Or decreasing the diameter would reduce the hydrocyclone capacity. Hence, a hydrocyclone design is optimum when all these factors, and their effects, produce a desirable balance. Gaining a deep understanding of the effect of the hydrocyclone geometrical and operational parameters on the flow field and its pertinent details is imperative to such optimum design. CFD has been promised to provide a reliable and affordable tool to do so. The effect of the above discussed design variables can be summarized as the followings. 1.3.7.1 Effect of Vortex Finder Diameter Decreasing vortex finder diameter leads to an increase in the residence time, and consequently results in a reduction in the cut size. Increasing vortex finder diameter increases the separation sharpness. A minimum of D/8 has been recommended (Bradley 1965) for the vortex finder 26 diameter, whereas a limit on a maximum for the vortex finder diameter recommends that it should not be greater than the LZVV so as not to disrupt the normal pattern of the inward radial flow (Bradley and Pulling, 1959). While the LZVV was mentioned to be at 0.43D, the overflow diameters in practical designs range from 0.16D to 0.5D (Svarovsky, 1984). 1.3.7.2 Effect of Cone Length and Angle While the cylindrical section can be very short in some hydrocyclones, or sometimes eliminated altogether, the conical section is imperative to the hydrocyclone functionality for a number of reasons. The momentum loss in the downward flow is balanced out with the reduction in the cross sectional area due to the narrowing passage in the cone. This narrowing is responsible for the similarity of the tangential velocity profiles across the hydrocyclone. Moreover, the cone diminishes secondary recirculation flows, thus enhancing the separation performance. Increasing the cone length increases the separation efficiency as well as the overall capacity. In practical designs the overall cyclone length varies between 2.7D to 7.7D, with small hydrocyclones taking the higher values. Decreasing the cone angle improves the separation efficiency by reducing the cut size. However, increasing the cone angle increases the separation sharpness. The conventional cono-cylindrical hydrocyclones designs are usually categorized as narrow angle designs (cone angle ranging from 6 to 25 degrees) and wide angle designs (cone angle ranging from 25 to 180 degrees). Narrow-angle designs are suitable for fine particle separation. A narrow passage dampens recirculating flows, increasing hydrocyclone efficiency for the separation of fine particles. 1.3.7.3 Effect of Inlet The inlet size, and its position in the case of multiple inlets, plays a significant role in the hydrocyclone flow field and separation performance. The inlet opening size determines the inlet velocity when a fixed flow rate is desired. Whereas, the inlet opening size determines the flow rate when a fixed inlet velocity is desired. Therefore, the inlet size (and numbers) determines the capacity, and to a large extent, the tangential velocity field pattern inside the hydrocyclone. Generally speaking, smaller inlet sizes improve the separation efficiency and reduce the cut size, 27 at the expense of capacity. Hence, the hydrocyclone design should offer a practical compromise between a reasonable separation performance and a satisfactory capacity. Kelsall (1953) reported an optimum inlet size of 0.08D, which was much smaller than the recommended optimum inlet size of 0.28D (Rietema 1961a, 1961b). Bradley (1965) suggested a variation of 0.14D to 0.17D for the inlet size as a good compromise. In practice, the inlet design can have an inlet size ranging from 0.13D to 0.29D (Svarovsky, 1984). 1.4 Turbulence Modelling CFD simulations of the flow field within hydrocyclones have been intensified with the advancement of computational facilities over the past couple of decades. Turbulence modelling has been shown to be crucial to the fidelity of CFD predictions. Highly swirling fluid flow inside the hydrocyclone with curved streamlines, flow reversal, vortex break down, and flow separation near the wall generates turbulence anisotropy. Therefore, basic turbulence models such as two-equation k-ε and its modified variants, due to their inherent turbulence isotropy assumption, fail to accurately predict the hydrocyclone flow field, as shown by other authors (Slack et al., 2000; Nowakowski et al., 2004; Delgadillo and Rajamani, 2005; Narasimha et al., 2007). Prior studies used turbulence models ranging from simple algebraic models to advanced, computationally expensive LES (Boysan et al., 1982; Ghodrat et al., 2014b; He et al., 1999; Hsieh and Rajamani, 1991; Saidi et al., 2012; Slack et al., 2000; Vieira et al., 2011). Streamline curvature diminishes the transport of turbulence in the forced-vortex region (direction of increase of angular momentum with radius) whereas the turbulence transport is enhanced or augmented in the free-vortex region (the direction of decrease of angular momentum with radius) (Launder et al., 1977). Furthermore, turbulence dispersion can have a profound effect on particle separation. Turbulence can increase separation efficiency by enhancing mixing. However, the significance of turbulence may diminish when particle size, i.e. length scale, increases (Olson and Kerekes 2017). Thus, turbulence plays a critical role in separation of smaller scale particles while it does not have a substantial influence of large scale particles. 28 An overview of the evolution of CFD simulations of the hydrocyclone flow field was discussed in review studies (Nowakowski et al., 2004; Narasimha et al., 2007). The following presents an account of the performance of various turbulence models, used over the decades, to predict the hydrocyclone flow pattern. 1.4.1 Performance of Turbulence Models There is consensus that conventional two-equation models such as standard k-ε are incapable of accurately capturing strong swirls and highly curved streamlines (Pericleous and Rhodes, 1986; Sevilla and Branion, 1997; Ma et al., 2000; Delgadillo and Rajamani, 2005). This incapability is due to the isotropy assumption that these models make for the eddy viscosity (see section 2.2). However, turbulence in the hydrocyclone flow field is categorically anisotropic due to flow reversal and the strong swirling nature of the flow. Yet, due to their mathematical simplicity and low computational costs, these models were not fully abandoned. Numerous attempts were made to introduce new forms of the k-ε model by adjusting its constants so as to bring numerical results closer to experimental data. Early attempts at numerically modelling the hydrocyclone flow field were made using simple algebraic turbulence models. Boysan et al. (1982) proposed a two-dimensional, algebraic stress model for turbulence to predict the velocity field and particle distributions. Their model was limited to gas cyclones and was too complicated and inefficient time-wise. Also, it did not take into account changes in the slurry viscosity at higher pulp concentrations. Pericleous and Rhodes (1986) developed a steady, axisymmetric mathematical PHOENICS code to measure the pulp velocity profile as well as particle distribution. They used a modified Prandtl mixing length model to account for the effect of swirl and particle presence on turbulence. The velocity predictions were in agreement with the experimental data of Kelsall (1952). However, the model was not tested for variation of inlet conditions or other geometries. Hsieh (1988) indicated that the gradients of the axial and tangential velocity components, in the radial direction, were of the same order magnitude. Thus, a modified Prandtl mixing length model was proposed with two different mixing length constants, one for each direction of the momentum equation, to take into account turbulence anisotropy. The axisymmetric model was developed 29 based on the physics of the fluid flow inside the hydrocyclone. The predicted results were successfully validated against LDV-measured data. Despite the excellent experimental validation of the model, its development to a hydrocyclone design tool was challenging (Hsieh and Rajamani, 1991). Dai et al. (1999) used a modified k-ε turbulence model in an axisymmetric study. They adjusted the eddy viscosity coefficient, the generation constant and the dissipation constant such that the calculated results would better match LDV measurements. They reported that the relative turbulence is strongest in the axial direction; thus the corresponding eddy viscosity was lowered to reflect that, as well as turbulence anisotropy. Moreover, the generation constant was increased, and the dissipation constant was decreased to make the model more suitable for the existing high velocity gradients and strong, anisotropic turbulence inside the hydrocyclone. They observed a two-peak tangential velocity pattern inside the cylindrical section with a second peak closer to the wall, which was attributed to existing strong swirl at the inlet. Dyakowski and Williams (1993) proposed a revised approach to overcome the two-equation model limitations that arose from the curved streamlines, rotational strains, and inaccurate assumptions of anisotropy and local equilibrium for turbulence. The standard k-ε equations were coupled with equations for the Reynolds stresses (Saffman, 1974). This resulted in a set of six partial differential equations along with six algebraic equations for the Reynolds stresses. Their approach resolved the standard k-ε model’s limitation of anisotropic turbulent viscosity, and successfully modelled the non-linear interaction between the mean vorticity and mean strain rate fields. However, the model had been proposed, and tested, only for small-diameter hydrocyclones, i.e. 10 - 44 mm. Using an axisymmetric model, Malhotra et al. (1994) incorporated into the standard k-ε, based on near-wall calculations and numerical optimization, a new formulation for the dissipation term to better predict the flow behavior. The overall performance of the proposed model was deemed reasonable, despite differences between the predicted results and experimental data. He et al. (1999) solved a modified k-ε for both axisymmetric and three-dimensional geometries. A curvature correction term was added to the dissipation equation in the standard k-ε model, based on the recommendation of Launderet al. (1977), to take into account the effect of the curvature of streamlines. Simplicity in numerical implementation of the correction term provided the model with an advantage over modified versions with additional equations. Nevertheless, the modified 30 k-ε model showed excellent agreement with the available experimental data. They reported that the axisymmetric model predicted the flow pattern reasonably accurate, even though the three-dimensional model produced the more accurate results. The three-dimensional calculations showed the non-axisymmetrical flow regions existed at the top of the cylindrical body; the flow approached axisymmetry farther away from the inlet in other regions and particularly in the cone. In terms of separation performance, the three-dimensional model was superior as the axisymmetric model under-predicted the particles cut-size; the error grew larger when the particles density was close to that of the carrying fluid, i.e. water. Thus, the axisymmetric model was inadequate for particle separation performance; nonetheless, its flow field predictions were adequately accurate. ReNormalization Group (RNG) is another variant of the k-ε model which was employed by several researchers (Griffiths and Boysan, 1996; Ma et al., 2000). Using the RNG model for their three-dimensional calculation, Ma et al. (2000) admitted that a more comprehensive model such as RSM was a more accurate and suitable choice for the hydrocyclone flow field. They indicated that the accuracy of the turbulence model was not independent of the cyclone geometry. This meant that while a specific turbulence model gave reasonably accurate predictions for a certain cyclone geometry, its accuracy could be impaired when applied to a different geometry. They concluded that given the dependency of the fidelity of turbulence modelling on the cyclone geometry, there was a need for development of a more universal turbulence model which could produce a reliably accurate performance for a wider range of geometries. In terms of particle separation, the study reported that the fate of particles is substantially influenced by their initial positions in the inlet; particles entering from the outer side of the inlet plane were more likely to be separated by the cyclone. Even though the modified variants of k-ε models improved the CFD predictions, they were largely limited to specific cyclone geometries. That, combined with computer advancements, led many authors to turn their attention to physically more realistic, and mathematically more elaborate, turbulence models. Averous and Fuentes (1997) used the RSM in axisymmetric hydrocyclone calculations. They pointed out that the relative turbulence was strongest in the radial direction, which was contrary to Dai et al. (1999). They suggested that the existing limitations of standard hydrocyclone models 31 was due to the empirically-obtained constants, upon which the models rely to produce ‘usable’ results. They concluded that these constants were for standardized experimental measurements, so it was possible that they would not work accurately for more complex geometries or flow fields. Moreover, the study found that the fate of particles was not influenced by their initial position in the inlet tube; classification was independent of the particle’s initial position in the inlet duct. This finding was in direct contradiction of that of Ma et al. (2000). Slack et al. (2000) conducted a three-dimensional CFD study for the hydrocyclone flow field modelling using both the RSM and LES. They found that despite the close agreement between the RSM results and the LDV measurements inside the cylindrical body, the RSM under-predicted the tangential velocity farther away from the inlet, and particularly in the cone. They noted that while the RSM predictions were within “engineering accuracy” for a wide variety of engineering problems they did not provide universality for bounded swirling flows. The study showed that the LES model produced a slightly superior performance in matching the experimental data as well as capturing time-dependent vortex precisions. However, it required a much finer mesh and longer computational time. On the other hand, the RSM simulations on a much coarser grid gave a satisfactory performance in investigating the time-averaged flow field inside a hydrocyclone in detail. Thus, the RSM was inexpensive compared to LES but not as accurate in capturing the fine details of the flow field. Cullivan et al. (2003b) used two forms of RSM, namely linear and quadratic pressure-strain (LPS and QPS) models for three-dimensional modelling of water flow inside a hydrocyclone. The pressure-strain term in the RSM is responsible for the redistribution of normal Reynolds stresses and the reduction of shear Reynolds stresses (this will be discussed in details in section 2.2.2). They reported significant difference between the predictions of the LPS and QPS models. The research indicated that the QPS model revealed details of the flow field structures which the LPS model failed to capture. The same authors, in a separate study (Cullivan et al., 2003a), stated that the QPS model led to an improved turbulence anisotropy predictions, and consequently to more accurate velocity and pressure fields, compared to the LPS model. They remarked that the QPS model should be viewed as a lower bound for pressure-strain modelling in CFD simulations of the 32 hydrocyclone flow field, and possibly a range of other industrial processes involving strong swirl and substantial recirculation. Under-prediction of tangential velocity by the RSM has also been reported in other research works. In a very interesting CFD study, Brennan (2006) reported that regardless of the grid refinement, the RSM persistently under-predicted the tangential velocity, which he attributed to a “problem” with the RSM for the hydrocyclone flow field. Furthermore, no significant difference was observed between the results of the LPS and QPS models. Narasimha et al. (2006) reported that once the air-core was established, the predictions from the LPS and QPS models were very similar, an observation which was in contradiction to that of Cullivan et al. (2003a). The study indicated the need for adjusting the model’s constants so as to improve the velocity predictions. The same authors (Brennan et al., 2009) reported that the RSM, compared to LES, gave more satisfactory results for water only simulations. However, their grids were coarse which could be the culprit for inferior LES performance. The authors admitted that what happened in practice might be different than their simulations results. It was concluded that the choice of turbulence needed further investigation. The RSM constants were proposed empirically by curve-fitting exercises on experimental data for simple shear flows with certain types of boundary conditions to represent a wide range of turbulence scales (Speziale et al., 1991). Like any other empirical correlations, the accurate applicability of the RSM constants is limited to the conditions from which they were proposed. In other words, the universally accepted default values for the model’s constants were not introduced to specifically represent the strong swirl, streamline curvature, and vortex break down present in the hydrocyclone flow field. Adjusting RSM constants to enhance the accuracy of CFD predictions seems a viable and desirable option in order to address this shortcoming. The newly adjusted coefficients could help overcome the model’s deficiency in the numerical calculations of the hydrocyclone flow field, and possibly bounded, swirling flows as a whole. In an attempt to improve the RSM performance for the hydrocyclone flow field by adjusting the model’s constants, Brennan (2006) reduced the return-to-isotropy constant, 𝐶1, and increased the fast pressure term, 𝐶2, in the LPS model (see section 2.2.2.2.1 for more details on these constants). He chose the linear model for recalibration since it had only two constants as opposed to the QPS’ seven constants. 33 Decreasing 𝐶1 caused little change, while increasing 𝐶2 much beyond its default value resulted in satisfactory improvements in the tangential velocity predictions. This improved performance came at the cost of considerable oscillations in the axial velocity field with time. Brennan (2006) suggested that adjusting the fast term constants in the QPS model would be a “useful next step” in developing a hydrocyclone-specific turbulence model. A second moment closure turbulence model, the Reynolds Stress Model (RSM), has been shown to be a minimum requirement for accurately capturing turbulence anisotropy in the hydrocyclone flow field (Murthy and Bhaskar, 2012; Cullivan et al., 2003a; Cullivan et al., 2004; Kuang et al., 2012). In spite of the good agreement between RSM results and experimental data in these studies, the standard RSM has been shown to lack generality when applied to hydrocyclones with different geometrical features (Slack et al., 2000). Moreover, under-prediction of the swirl velocity has been reported by several authors (Slack et al., 2000; Brennan, 2006). The shortcomings of the RSM in capturing the details of the hydrocyclone flow field are due to the fact that the model’s empirical constants were obtained by curve fitting over a set of experimental data to represent a wide range of turbulence flow fields. Hence, the constants do not represent the bounded, highly swirling hydrocyclone flow field. It has been well established in the literature that amongst the available turbulence techniques in CFD, excluding DNS (direct numerical simulation), the LES turbulence model provides the most accurate details of the hydrocyclone flow field including the velocity profiles, turbulence fluctuations and intensity, pressure drop, and particle separation (Slack et al., 2000; Delgadillo and Rajamani, 2005; Brennan, 2006; Narasimha et al., 2006; Delgadillo and Rajamani, 2009). However, the LES requirements for a very fine grid combined with small time step size, makes the model highly expensive to employ, and acts as a deterrent for researchers and engineers to adopt it as a design tool. Therefore, there exists a need for developing an accurate and affordable turbulence model for the hydrocyclone flow field and for bounded swirling flows in general. A modified version of the RSM could fulfill the demand and provide a great compromise between accuracy and affordability. 34 1.4.2 Effect of Swirl on Turbulence In an attempt to better understand the hydrocyclone flow field, we present a summary of the established literature for a more general phenomenon in fluid dynamics: the interaction between swirl and turbulence. This interaction occurs in many industrial applications including removal or classification of particles by cyclone separators in the pulp and paper, food, chemical and mineral processing industries (Nowakowski and Dyakowski, 2003; He et al., 1999). The behavior of turbulence in swirling flow fields has been attractive to researchers for the past several decades due to its engineering applications such as in cyclone separators, turbomachinery, heat exchangers, and turbulence modeling among others (Jones and Launder, 1972; Howard et al., 1980; Murakami and Kikuyama, 1980). A rotating pipe presents a classic fluid dynamics problem that can be very useful for improving turbulence models for swirling flows, including the hydrocyclone flow field (Orlandi and Fatica, 1997). It has been established by numerous studies that pipe rotation dampens the turbulence leading to the re-emergence of laminar characteristics, also referred to as ‘laminarization’, ‘reverse-transition’, ‘inverse-transition’, or ‘re-laminarization’. One of the early attempts to understand the effect of swirl on turbulence was made by White (1964). In an experimental study on the subject of fluid flow inside a rotating pipe, he reported opposite responses of the flow to the pipe rotation depending on the flow regime. Rotation in the laminar regime led to an increase in the pressure loss. On the contrary, the pipe rotation caused a considerable reduction in the pressure drop in the turbulent regime. However, the pressure loss remained rather unaffected by rotation, for flows with Reynolds numbers higher than 8000. Moreover, suppression of turbulence momentum exchange due to the pipe rotation was experimentally observed. The author admitted that, despite a dimensionless representation, the results had been obtained for only one pipe size; thus, the possible effects of scaling had not been considered. Murakami and Kikuyama (1980) conducted an experimental investigation on the flow pattern and hydraulic loss in a rotating pipe, with the pipe length varying from 30D to 160D, where D is the pipe diameter. They reported that the degree to which turbulence was suppressed was a function of the axial distance from the inlet and the rotation number, 𝑁 = 𝑅𝜔 𝑈𝑚⁄ (inverse of Rossby 35 number used in geophysical problems), defined as the ratio of the rotating wall velocity (R is the pipe radius and ω is the pipe rotational speed) to the flow mean axial velocity (Um). They surprisingly reported that at a high rotation number of N=8.6 the flow showed neither an entirely laminar axial profile nor a wholly forced-vortex tangential component distribution. Opposing effects of the rotation on the turbulence in different regions of the pipe was also reported in other studies (Kikuyama et al., 1983; Nishibori et al., 1987). Using a rotated hot wire method, Kikuyama et al. (1983) measured the mean velocity components and pressure drop. They reported opposing effects of rotation on turbulence in different regions of the pipe. The ‘intensive’ shear stress caused by the wall rotation increased the shear Reynolds stresses, which in turn led to a destabilizing effect in the vicinity of the inlet, within the entrance length (x/D < 10). Farther away from the inlet (x/D > 10), however, the swirling flow developed, and the resulting centrifugal force suppressed the turbulence leading to a stabilizing effect. The study concluded that the turbulence intensity would eventually drop below that in the stationary pipe condition. Nishibori et al. (1987) conducted an experimental study on the laminarization in a rotating pipe for rotation numbers of N = 0, 1, and 3. They reported a ‘complicated’ and alternatingly changing flow pattern in the inlet region. At all N’s the axial velocity distribution was flattened near the pipe centre of the inlet regions. They attributed the flatness to the turbulence suppression by the centrifugal force of the swirling flow, which transformed the axial velocity profile in the boundary layer to a laminar-like boundary layer velocity profile. Further downstream, however, the axial velocity profile started to deform into a parabola and gradually took the laminar shape. They observed that, with visualization techniques using aluminum flakes, a destabilizing effect caused by the pipe rotation enhanced turbulence in the laminarized boundary layer and impeded the laminarization. The stabilizing effect of the centrifugal forces, however, grew dominant at sections further downstream of the inlet for the same rotation number, or with an increase of the rotation number (to N=3) at the same section. Additionally, they found that laminarization is enhanced, at the same rotation number, at the inlet region by decreasing the axial Reynolds number. Reich and Beer (1989) experimentally and analytically investigated the fluid flow and heat transfer in a rotating pipe for a range of Reynolds and rotation numbers. They also reported a considerable reduction in the flow resistance and heat transfer in the pipe due to rotation. Their velocity and 36 temperature profiles indicated a marked laminarization. The suppression of turbulent migration of fluid particles, due to the developed centrifugal force, was responsible for the observed laminarization and drag reduction. Their observations and conclusions were consistent with those of Nishibori et al. (1987). Orlandi and Fatica (1997) employed DNS to study the flow structure in a rotating pipe with the focus being on predicting and explaining the drag reduction phenomenon in a rotating pipe, and on the near-wall turbulence. The highest rotation number used in the study was N =2 and the pipe length was L= 7.5D. The authors reported similar findings as previous studies; rotation resulted in drag reduction and at high rotation numbers the mean axial velocity profile took the parabolic laminar form. They claimed that the drag reduction occurred due to modifications of the near-wall vertical structure. The effects of the modifications could either be in addition to the centrifugal force effect reported before (Nishibori et al., 1987; Reich and Beer, 1989), or the effect of centrifugal force manifested in the modifications of the vertical structure. Feiz et al. (2003) used the LES model for a rotating pipe with L=10D, and two rotation numbers of N = 1 and 2. The authors reported a decrease in turbulence intensity with an increase in the pipe rotational speed. This observation was attributed to the stabilizing effect of the centrifugal force. Moreover, it was reported that the dynamic model gave superior performance over the Smagorinsky model. The difference between these two LES models will be discussed in detail in section 2.2.3. In summary, the effect of rotation on turbulence has been investigated intensively via experimental, numerical and analytical investigations over the years. The effect of the centrifugal forces on drag reduction, turbulence intensity, and Reynolds stresses, has been closely studied. There is consensus in the literature that the pipe rotation generates a stabilizing effect downstream of the inlet, which suppresses turbulence and leads to the laminarization phenomena. In regards to turbulence modelling, LES and DNS are the most advanced and accurate numerical models to tackle these types of swirling turbulence problems. However, they incur very high computational costs and are not affordable to be used for industrial applications. The RSM has been shown to be a reasonable compromise between accuracy and cost. Specifically, it has been shown to give very accurate predictions for the velocity field and turbulence in a rotating pipe (Reich and Beer, 1989). 37 Thus, RSM remains a viable design tool for engineering applications owing to the high computational costs incurred by the more accurate LES and DNS techniques. 1.5 Summary of Literature and Objectives Numerous studies have been conducted to understand the complex flow patterns of swirling flows in general, and inside hydrocyclones in particular. Several authors have investigated the interaction between swirl and turbulence. Many other authors have performed studies to understand the hydrocyclone flow field. However, there are no studies, to the best knowledge of the author, connecting the two problems of the very same class of turbulent flows. It is beneficial for increasing our understanding of the hydrocyclone flow field pattern to bridge the two investigated fields and initiate an interdisciplinary research approach. Such an approach will be useful in developing a more accurate turbulence model specifically for swirling flows and hydrocyclones. There is consensus in the literature that there is a need to develop a new turbulence model which provides generality for the flow field inside a hydrocyclone, irrespective of its geometry. One way to achieve such a model is to fine-tune the correlated constants of the more accurate turbulence models that exist. Attempts have been made to do so, with reasonable success, for the linear pressure-strain RSM. However, the quadratic pressure-strain model has been reported by several researchers to produce superior performance than the linear model. Nevertheless, the high computational costs of fine-tuning the quadratic model, due to considerably fewer constants of the LPS, has proved to be a deterrent in attempts to customize the QPS model for the hydrocyclone flow field. In this work, a sensitive study is performed on the standard constants of the QPS model with the objective of finding optimally accurate numerical values for the constants. The effect of the hydrocyclone design and flow variables on its separation performance has received a great deal of attention over the years. Much less attention has been devoted to the influence of those variables on the hydrodynamics of the flow other than the velocity field. Generally speaking, a vast majority of studies limited their flow field analysis to examine the effect of design variables on the velocity field. In other words, they stopped short of a more in-depth analysis of the interplay between different forces in the hydrocyclone flow field. They also failed to display how that interplay is affected by changes in the design variables, and subsequently how 38 these changes in hydrodynamics could be linked with changes in separation performance. In this study, a closer look is taken at the significance of different forces in the momentum balance, and their dependency on the design and flow variables. The objectives of this research work is summarized as follows, 1. There is a critical rotation number beyond which the convective forces become the major competitor for the pressure force, outgrowing the turbulence forces in the process. If we can establish a reliable threshold for that critical rotation number, provided that such number is universal regardless of other flow conditions, general guidelines could be provided for modelling turbulence in swirling flows. The hope would be that finding such criterion, and consequently producing guidelines, would enable engineers to employ simpler models than the RSM or LES for swirling flows with rotation numbers lower than the would-be-found critical rotation number. 2. The hydrodynamics of the hydrocyclone flow field is immensely dependent on the design parameters. By performing a force balance analysis, as opposed to a kinematic analysis, useful dimensionless groups could be proposed for scaling hydrocyclones, and predicting their flow field, force balance interplay, and ultimately separation performance. 3. The standard quadratic pressure-strain RSM is a reasonable and satisfactory compromise between accuracy and computational costs for the hydrocyclone flow field. However, it was not specifically designed for predicting swirling flows with curved streamlines. Some of its universal constants might be more influential than others in numerical predictions of the hydrocyclone flow field. Correlating those constants to more accurately represent the swirl and streamline curvature existing in the hydrocyclone will be an enormously useful step towards developing an affordable, yet accurate, hydrocyclone-specific turbulence model. Such constants can be identified by a sensitive study (changing one constant at a time) to understand the extent of their influence on the model’s predictions. Constants can be categorized into strong influence, moderate influence, and zero influence groups. Identifying strong influence, and to a lesser extend medium influence, constants will help facilitate designing a hydrocyclone-specific model. Consequently, we will be able to fine-tune the constants accordingly for optimal accuracy and cost-effectiveness. 39 Chapter 2: Methodology 2.1 Governing Equations The Navier-Stokes equations of conservation of mass (continuity) and transport of linear momentum for incompressible, Newtonian viscous flows with constant properties, are the governing equations of the flow fields studied in this text. These equations are expressed, in their conservative form, in the Cartesian coordinate system as follows, Continuity 𝜕𝑢𝑖𝜕𝑥𝑖= ∇⃗ ∙ ?⃗? = 0 (2.1) Momentum 𝜌𝜕𝑢𝑖𝜕𝑡+ 𝜌𝜕𝜕𝑥𝑗(𝑢𝑗𝑢𝑖) = −𝜕𝑝𝜕𝑥𝑖+𝜕𝜏𝑖𝑗𝜕𝑥𝑗+ 𝜌𝑔𝑖 (2.2) 𝑢𝑖 is the instantaneous velocity vector, 𝑥𝑖 is the position, ρ is density, t is time, p is instantaneous pressure, 𝑔𝑖 is the gravitational body force, and 𝜏𝑖𝑗 is the instantaneous stress tensor, defined as, 𝜏𝑖𝑗 = 2𝜇𝑠𝑖𝑗 (2.3) where μ is the fluid molecular viscosity and 𝑠𝑖𝑗 is the strain-rate tensor, defined for an incompressible Newtonian fluid with constant properties by, 𝑠𝑖𝑗 =12(𝜕𝑢𝑖𝜕𝑥𝑗+𝜕𝑢𝑗𝜕𝑥𝑖) (2.4) In the Cartesian coordinate the stress tensor takes the following form, 40 𝜏𝑖𝑗 = [𝜏𝑥𝑥 𝜏𝑥𝑦 𝜏𝑥𝑧𝜏𝑦𝑥 𝜏𝑦𝑦 𝜏𝑦𝑧𝜏𝑧𝑥 𝜏𝑧𝑦 𝜏𝑧𝑧] =( 2𝜇𝜕𝑢𝑥𝜕𝑥𝜇 (𝜕𝑢𝑥𝜕𝑦+𝜕𝑢𝑦𝜕𝑥) 𝜇 (𝜕𝑢𝑥𝜕𝑧+𝜕𝑢𝑤𝜕𝑥)𝜇 (𝜕𝑢𝑦𝜕𝑥+𝜕𝑢𝑥𝜕𝑦) 2𝜇𝜕𝑢𝑦𝜕𝑦𝜇 (𝜕𝑢𝑦𝜕𝑧+𝜕𝑢𝑧𝜕𝑦)𝜇 (𝜕𝑢𝑤𝜕𝑥+𝜕𝑢𝑥𝜕𝑧) 𝜇 (𝜕𝑢𝑧𝜕𝑦+𝜕𝑢𝑦𝜕𝑧) 2𝜇𝜕𝑢𝑧𝜕𝑧 ) (2.5) Note that all the above-mentioned variables take their instantaneous magnitudes. Now, we proceed by decomposing (Reynolds decomposition) the turbulent quantity into the sum of the mean value (denoted in uppercase) of the quantity and its fluctuating value (denoted by the prime sign). Hence, 𝑢𝑖 = 𝑈𝑖 + 𝑢𝑖′, 𝑝 = 𝑃 + 𝑝′, 𝑠𝑖𝑗 = S𝑖𝑗 + 𝑠𝑖𝑗′ , and 𝜏𝑖𝑗 = Τ𝑖𝑗 + 𝜏𝑖𝑗′ . Substituting these quantities into equations (2.1) - (2.5) and taking ensemble (time) average of the instantaneous Navier-Stokes equations, i.e. equations (2.1) and (2.2), results in 𝜌𝜕𝑈𝑖𝜕𝑡+ 𝜌𝑈𝑗𝜕𝑈𝑖𝜕𝑥𝑗= −𝜕𝑃𝜕𝑥𝑖+𝜕𝑇𝑖𝑗𝜕𝑥𝑗+𝜕𝜏𝑖𝑗𝑡𝜕𝑥𝑗+ 𝜌𝑔𝑖 (2.6) where 𝑇𝑖𝑗 is modeled from the viscous fluid model and 𝜏𝑖𝑗𝑡 from the turbulence model. They are defined as follows, 𝑇𝑖𝑗 = 2𝜇𝑆𝑖𝑗 = 𝜇 (𝜕𝑈𝑖𝜕𝑥𝑗+𝜕𝑈𝑗𝜕𝑥𝑖) (2.7) and 𝜏𝑖𝑗𝑡 = −𝜌𝑢𝑖′𝑢𝑗′̅̅ ̅̅ ̅̅ (2.8) Equation (2.6) is known as the Reynolds-Averaged Navier-Stokes (RANS) equation. The quantity −𝜌𝑢𝑖′𝑢𝑗′̅̅ ̅̅ ̅̅ is referred to as the Reynolds stresses and 𝜏𝑖𝑗𝑡 is known as the Reynolds stress tensor. Substituting equation (2.7) into equation (2.6) and using the continuity equation leads to the RANS equations for fluid flows with constant properties, 41 𝜌𝜕𝑈𝑖𝜕𝑡+ 𝜌𝑈𝑗𝜕𝑈𝑖𝜕𝑥𝑗= −𝜕𝑃𝜕𝑥𝑖+ 𝜇𝜕2𝑈𝑖𝜕𝑥𝑗2 +𝜕𝜏𝑖𝑗𝑡𝜕𝑥𝑗+ 𝜌𝑔𝑖 (2.9) Calculating the Reynolds stress tensor presents the “fundamental problem” of turbulence: “in order to compute all mean-flow properties of the turbulent flow under consideration, we need a prescription for computing 𝑢𝑖′𝑢𝑗′̅̅ ̅̅ ̅̅ ”(Wilcox, 2006). Note that the Reynolds stress tensor is symmetric (𝜏𝑖𝑗𝑡 = 𝜏𝑗𝑖𝑡 ). Thus, the Reynolds averaging process resulted in six additional unknowns and no new equations. Therefore, in addition to the mean pressure and the three mean velocity components for three-dimensional flows, we now have six unknown Reynolds stresses. This presents a set of four equations (one continuity and three momentum equations) with ten unknowns. Hence, the system of equations is not closed; we must introduce enough equations to be able to solve for our unknowns. The problem of mathematically closing the turbulence model is usually referred to as the turbulence closure problem. In order to introduce additional equations, we will derive equations for 𝑢𝑖′𝑢𝑗′ by taking moments of the Navier-Stokes equations; multiply equation (2.2) for 𝑢𝑖 and 𝑢𝑗 , by 𝑢𝑗′ and 𝑢𝑖′, respectively, sum up the results of the two operations, and take the ensemble average of the total summation. Following this procedure, combined with some re-arrangement and manipulations, leads to the Reynolds-stress equations defined as, 𝜕𝜏𝑖𝑗𝑡𝜕𝑡+ 𝐶𝑖𝑗 = 𝐷𝑇,𝑖𝑗 +𝐷𝐿,𝑖𝑗 + 𝑃𝑖𝑗 − ∅𝑖𝑗 + 𝜀𝑖𝑗 (2.10) 𝐶𝑖𝑗 =𝜕𝜕𝑥𝑘(𝑈𝑘𝜏𝑖𝑗𝑡 ) (2.11) 𝐷𝑇,𝑖𝑗 =𝜕𝜕𝑥𝑘(𝜌𝑢𝑖′𝑢𝑗′𝑢𝑘′ + 𝑝′(𝛿𝑘𝑗𝑢𝑖′ + 𝛿𝑖𝑘𝑢𝑗′)) (2.12) 𝐷𝐿,𝑖𝑗 = −𝜕𝜕𝑥𝑘(𝜈𝜕𝜏𝑖𝑗𝑡𝜕𝑥𝑘) (2.13) 42 𝑃𝑖𝑗 = (𝜏𝑖𝑘𝑡𝜕𝑈𝑗𝜕𝑥𝑘+ 𝜏𝑗𝑘𝑡 𝜕𝑈𝑖𝜕𝑥𝑘) (2.14) Φ𝑖𝑗 = 𝑝′ (𝜕𝑢𝑖′𝜕𝑥𝑗+𝜕𝑢𝑗′𝜕𝑥𝑖) (2.15) 𝜀𝑖𝑗 = 2𝜇𝜕𝑢𝑖′𝜕𝑥𝑘𝜕𝑢𝑗′𝜕𝑥𝑘 (2.16) The physical descriptions of the above terms of the Reynolds stress equation are as follows, 𝐶𝑖𝑗 is the transport of the momentum (Reynolds stresses) by convection. 𝐷𝑇,𝑖𝑗 represents the transport of the momentum (Reynolds stresses) by turbulent diffusion. 𝐷𝐿,𝑖𝑗 represents the transport of the momentum (Reynolds stresses) by molecular diffusion. 𝑃𝑖𝑗 represents the production of the stresses through the interaction of the Reynolds stresses with the gradient of the mean flow velocity field. Φ𝑖𝑗 is called the pressure-strain (PS) term and represents the tendency of the flow to restore isotropy. 𝜀𝑖𝑗 is viscous dissipation of the Reynolds stresses. The turbulence closure problem is solved by devising correlations, based on the physics and in terms of known flow properties, such that new equations are introduced to approximate the unknowns. In doing so, we will close the system of equations. This is, in essence, the function of turbulence modelling. 2.2 Turbulence Modeling The following presents details and assumptions of the turbulence models used in this study, namely the standard k-ε model, the Reynolds Stress Model, and the Large Eddy Simulation model. 43 2.2.1 Standard k-ε Model The k-ε model is the most widely used and validated turbulence model in engineering. It was developed based on the fundamental presumption that the action of the Reynolds stresses on the mean flow are analogous to those of the viscous stresses. The model is categorized as a turbulence energy equation model. The kinetic energy (per unit mass) of the turbulent fluctuations, referred to as the turbulence kinetic energy, k, is the basis for this model. It represents the specific turbulence kinetic energy, and is defined by, 𝑘 =12𝑢𝑖′𝑢𝑖′̅̅ ̅̅ ̅̅ =12(𝑢1′2̅̅ ̅̅ + 𝑢2′2̅̅ ̅̅ + 𝑢3′2̅̅ ̅̅ ) (2.17) Note that there exists a proportionality between the trace of the Reynolds stress tensor and the turbulent kinetic energy, 𝜏𝑖𝑖𝑡 = −𝜌𝑢𝑖′𝑢𝑖′̅̅ ̅̅ ̅̅ = −2𝜌𝑘 (2.18) This proportionality helps us derive an equation for k by taking the trace of the Reynolds stress equation; use equation (2.10) for i=j=1, i=j=2, and i=j=3, and sum the three equations. This procedure results in the following transport equation for the turbulence kinetic energy, 𝜌𝜕𝑘𝜕𝑡+ 𝜌𝜕𝜕𝑥𝑗(𝑈𝑗𝑘) = 𝜏𝑖𝑗𝑡 𝜕𝑈𝑖𝜕𝑥𝑗+𝜕𝜕𝑥𝑗(𝜇𝜕𝑘𝜕𝑥𝑗−12𝜌𝑢𝑖′𝑢𝑖′𝑢𝑗′̅̅ ̅̅ ̅̅ ̅̅ − 𝑝′𝑢𝑗′̅̅ ̅̅ ̅̅ ) − 𝜀 (2.19) where 𝜀 is the turbulence dissipation rate (per unit mass), and is defined as, 𝜀 = 𝜇𝜕𝑢𝑖′𝜕𝑥𝑘𝜕𝑢𝑖′𝜕𝑥𝑘̅̅ ̅̅ ̅̅ ̅̅ ̅̅ (2.20) It is worth mentioning that the trace of the pressure-strain tensor of equation (2.10), Φ𝑖𝑗, vanishes for incompressible flows (through using the continuity equation). The physics behind different terms in equation (2.19) is explained below. The unsteady and convection terms on the left-hand side (LHS) are the rate of change of k and transport of k by convection, respectively. The first term on the right-hand side (RHS) is called the 44 production term, and defines the rate of production of k. In other words, it represents the rate at which the mean flow transfers kinetic energy to the turbulence. Equivalently, the production term can be seen as the rate of work done by the mean strain rate against the turbulent stresses. The second term on the RHS, i.e. 𝜇𝜕2𝑘𝜕𝑥𝑗2, is the known as molecular diffusion, and represents the transport of k due to the molecular diffusion. The third term on the RHS, i.e. the triple velocity term, is referred to as turbulent transport, and represents transport of k by turbulent fluctuations. The fourth term on the RHS represents transport of k by correlation of turbulent fluctuations of pressure and velocity, and is called pressure diffusion. Finally, the dissipation term represents the rate at which the fluctuating strain rate does work against the fluctuating viscous stresses; it is equivalent to the rate at which the turbulence energy is lost through conversion into thermal energy, increasing the fluid internal energy. Mathematically speaking, the unsteady, convection, and molecular diffusion terms are exact, therefore they do not require modelling. However, the production, turbulent transport, pressure diffusion, and dissipation terms involve unknown correlations. Thus, 𝜏𝑖𝑗𝑡 , the turbulent transport, pressure diffusion, and dissipation, must be defined by relevant approximations or models to close the set of equations. 2.2.1.1 Reynolds Stress Tensor The k-ε model employs Boussinesq approximation to model the Reynolds-stress tensor, which assumes the Reynolds stresses are proportional to the mean deformation rates, and is defined as, 𝜏𝑖𝑗𝑡 = 2𝜇𝑡𝑆𝑖𝑗 −23𝜌𝑘𝛿𝑖𝑗 (2.21) where, as mentioned in equation (2.7), S is the mean strain-rate tensor, 𝛿𝑖𝑗is the Kronecker delta (𝛿𝑖𝑗 = 1 if 𝑖 = 𝑗, and 𝛿𝑖𝑗 = 0 when 𝑖 ≠ 𝑗), and 𝜇𝑡 is turbulent eddy viscosity, defined as, 𝜇𝑡 = 𝐶𝜇𝜌𝑘2𝜀 (2.22) with the closure coefficient 𝐶𝜇 = 0.09. 45 2.2.1.2 Turbulent Transport and Pressure Diffusion In analogy to molecular transport processes, the gradient-diffusion approximation is made to account for turbulent transport of scalar quantities, e.g. −𝑢𝑗′𝜑′̅̅ ̅̅ ̅̅ ≈ 𝜈𝑇 ∂Ψ/ ∂𝑥𝑗. However, there does not exist such analogy for the pressure diffusion term. With no reliable experimental data (Wilcox, 2006), it was assumed that the pressure-diffusion term and the turbulent transport combined, behave like a gradient-transport phenomenon. Hence, the sum is assumed to be, 12𝜌𝑢𝑖′𝑢𝑖′𝑢𝑗′̅̅ ̅̅ ̅̅ ̅̅ + 𝑝′𝑢𝑗′̅̅ ̅̅ ̅̅ = −𝜇𝑡𝜎𝑘𝜕𝑘𝜕𝑥𝑗 (2.23) where 𝜎𝑘 is a closure coefficient and is equal to unity. Substituting equation (2.23) into equation (2.19) leads to, 𝜌𝜕𝑘𝜕𝑡+ 𝜌𝜕𝜕𝑥𝑗(𝑈𝑗𝑘) = 𝜏𝑖𝑗𝑡 𝜕𝑈𝑖𝜕𝑥𝑗+𝜕𝜕𝑥𝑗[(𝜇 +𝜇𝑡𝜎𝑘)𝜕𝑘𝜕𝑥𝑗] − 𝜀 (2.24) Equation (2.24) is the modelled form of the transport of turbulent kinetic energy, i.e. equation (2.19), and is used in all variants of the k-ε model including the standard k-ε. In order to close the set of the equations we still need to define proper correlations for the dissipation term, ε. The exact equation for the dissipation rate is obtained by taking the following three steps: first take derivative of equation (2.2) for 𝑢𝑖, with respect to 𝑥𝑗, i.e. 𝜕/𝜕𝑥𝑗 , then multiply the derivative by 𝜕𝑢𝑖′/𝜕𝑥𝑗, and finally take time average of the resulting equation and multiply the time-averaged equation by 2𝜈. Performing a substantial amount of algebraic manipulation on the final time-averaged equation, leads to the exact equation for ε (Wilcox, 2006), 46 𝜕𝜀𝜕𝑡+𝜕𝜕𝑥𝑗(𝑈𝑗𝜀)= −2𝜈(𝑢𝑖,𝑘′ 𝑢𝑗,𝑘′̅̅ ̅̅ ̅̅ ̅̅ ̅ + 𝑢𝑘,𝑖′ 𝑢𝑘,𝑗′̅̅ ̅̅ ̅̅ ̅̅ ̅)𝜕𝑈𝑖𝜕𝑥𝑗− 2𝜈 𝑢𝑘′ 𝑢𝑖,𝑗′̅̅ ̅̅ ̅̅ ̅𝜕2𝑈𝑖𝜕𝑥𝑘𝜕𝑥𝑗− 2𝜈 𝑢𝑖,𝑘′ 𝑢𝑖,𝑚′ 𝑢𝑘,𝑚′̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅− 2𝜈2𝑢𝑖,𝑘𝑚′ 𝑢𝑖,𝑘𝑚′̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅ +𝜕𝜕𝑥𝑗(𝜈𝜕𝜀𝜕𝑥𝑗− 𝜈𝑢𝑗′𝑢𝑖,𝑚′ 𝑢𝑖,𝑚′̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅ − 2𝜈𝜌𝑝,𝑚′ 𝑢𝑗,𝑚′̅̅ ̅̅ ̅̅ ̅̅ ̅) (2.25) It is evident that the exact equation for ε is substantially more complex than that for k, i.e. equation (2.19). The task of finding proper closure approximations seems to be impossible to complete by experimentalists due to the highly complex nature of the equation. Despite some insightful DNS studies (Mansour et al., 1988) for the low-Reynolds-number flows, there is little hope of establishing suitable approximations to close the exact equation for the dissipation rate (Wilcox, 2006). After some physical considerations, the modelled version equation for the dissipation rate assumes the following form, 𝜌𝜕𝜀𝜕𝑡+ 𝜌𝜕𝜕𝑥𝑗(𝑈𝑗𝜀) = 𝐶𝜀1𝜀𝑘𝜏𝑖𝑗𝑡 𝜕𝑈𝑖𝜕𝑥𝑗− 𝐶𝜀2𝜀2𝑘+𝜕𝜕𝑥𝑗[(𝜇 +𝜇𝑡𝜎𝜀)𝜕𝜀𝜕𝑥𝑗] (2.26) where the closure coefficients take the following numerical values: 𝐶𝜀1 = 1.44, 𝐶𝜀2 = 1.92, and 𝜎𝜀 = 1.3. The standard k-ε model (Launder and Sharma, 1974) is formed by the modelled equations for the eddy viscosity, turbulent kinetic energy, and turbulent dissipation rate, i.e. equations (2.22), (2.24), and (2.26), along with their specified closure coefficients. 2.2.2 Reynolds Stress Model The Reynolds Stress Model solves for the Reynolds stresses transport equations along with an equation for the dissipation rate. Unlike the k-ε (or k-ω) model, the RSM does not make the isotropic eddy-viscosity assumption, i.e. it goes beyond Boussinesq approximation by solving transport equations for the Reynolds stresses, rather than simplifying them. The Reynolds stress equation, i.e. equation (2.10), provides the exact differential equation for the transport of the Reynolds stress, and lays the foundation for the Reynolds Stress Model. The convection (𝐶𝑖𝑗), 47 molecular diffusion (𝐷𝐿,𝑖𝑗), and production terms (𝑃𝑖𝑗), are exact and do not require modelling; they are retained in the Reynolds stress transport equation of (2.10), as stated in equations (2.11), (2.13), and (2.14). However, the turbulent diffusion, dissipation rate, and pressure-strain terms require modelling in order to close the system of equations. The modellings favoured by most commercial CFD codes, including ANSYS Fluent, are explained in the following subsections. Comprehensive details about various models for these terms have been given elsewhere (Launder et al., 1975; Rodi, 1980; Chen and Jaw, 1998; Wilcox, 2006). 2.2.2.1 Turbulent Diffusion Transport The most widely used method in modelling 𝐷𝑇,𝑖𝑗 is to assume a gradient-diffusion (gradient-transport) process. A simple approximation was proposed by Daly and Harlow (1970), 𝐷𝑇,𝑖𝑗 = 𝐶𝑠𝜕𝜕𝑥𝑘(𝑘𝜏𝑘𝑙𝑡𝜌𝜀𝜕𝜏𝑖𝑗𝑡𝜕𝑥𝑙) (2.27) where 𝐶𝑠 is a closure coefficient. Despite its mathematical simplicity, the approximation does not produce a symmetric tensor in 𝐷𝑇,𝑖𝑗, therefore rendering it inconsistent with the symmetrical nature of the tensor in its indices (Wilcox, 2006). Moreover, the approximation can result in numerical instability (Ansys, 2011); hence a simplified scalar version of the turbulent diffusion transport (Lien and Leschziner, 1994) is used in which the rate of transport of Reynolds stresses by diffusion and the gradient of Reynolds stresses are assumed to be proportional (Versteeg and Malalasekera, 2007), 𝐷𝑇,𝑖𝑗 =𝜕𝜕𝑥𝑘(−𝜈𝑡𝜎𝑘𝜕𝜏𝑖𝑗𝑡𝜕𝑥𝑘) (2.28) where 𝜈𝑡 = 𝜇𝑡 𝜌⁄ is the turbulent kinematic eddy viscosity (𝜇𝑡 defined by equation (2.22) ) and 𝜎𝑘=0.82 (Lien and Leschziner, 1994), hence different than 𝜎𝑘=1 in the standard k-ε model, i.e. equation (2.23). 48 2.2.2.2 Pressure-Strain Term The pressure-strain (PS) term plays a crucial role in the accuracy of the RSM predictions and its appropriate modelling is crucial to the fidelity of CFD simulations. Moreover, owing to the fact that unmeasurable correlations exist in Φ𝑖𝑗, it demands enormous physical insight into the problem to propose a reasonable and suitable closure approximation. Thus, the pressure-strain term has received the greatest deal of attention from turbulence researchers. To calculate the pressure fluctuations, 𝑝′, the following equation should be solved, 1𝜌∇2𝑝′ = − [𝜕2𝜕𝑥𝑖𝜕𝑥𝑗(𝑢𝑖′𝑢𝑗′ − 𝑢𝑖′𝑢𝑗′̅̅ ̅̅ ̅̅ ) + 2𝜕𝑈𝑖𝜕𝑥𝑗𝜕𝑢𝑗′𝜕𝑥𝑖] (2.29) This Poisson’s equation for the pressure fluctuations is obtained by time averaging the difference between the time-averaged Navier-Stokes equations and the instantaneous Navies-Stokes equations (first subtract the time-averaged Navier-Stokes from the instantaneous one, then time average). There are two general approaches to the pressure-strain modelling problem, usually referred to as the linear pressure-strain (LPS) and the quadratic pressure-strain (QPS) models. 2.2.2.2.1 Linear Pressure-Strain Integrating equation (2.29) over a large volume (vol), then applying Green’s theorem, and finally multiplying by (𝜕𝑢𝑖′𝜕𝑥𝑗+𝜕𝑢𝑗′𝜕𝑥𝑖), before time averaging results in the following expression for the PS, (PS) = Φ𝑖𝑗 = 𝑝′ (𝜕𝑢𝑖′𝜕𝑥𝑗+𝜕𝑢𝑗′𝜕𝑥𝑖)=14𝜋∫ 𝜌 { [𝜕2(𝑢𝑘′ 𝑢𝑙′)𝜕𝑥𝑙𝑥𝑚]∗(𝜕𝑢𝑖′𝜕𝑥𝑗+𝜕𝑢𝑗′𝜕𝑥𝑖)̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ vol+ 2(𝜕𝑈𝑙𝜕𝑥𝑚)∗(𝜕𝑢𝑚′𝜕𝑥𝑙)∗(𝜕𝑢𝑖′𝜕𝑥𝑗+𝜕𝑢𝑗′𝜕𝑥𝑖)̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ } 𝑑 vol|𝐱 − 𝐲| (2.30) 49 where the integration is carried out over the entire y space and values with a star in the integrand are at y, and values without a star in the integrand are at x. For more details on the derivation of equation (2.30) refer to (Launder et al., 1975; Chen and Jaw, 1998; Wilcox, 2006). Note that the first term of the integrand (in the bracket) of equation (2.30) explicitly involves only the fluctuating velocity components, whereas the second term involves the mean velocity gradient (mean strain rate). Thus, the first term is usually referred to as the slow PS or the return-to-isotropy term and is denoted by Φ𝑖𝑗,1, and the second term is known as the rapid (or fast) PS and is shown by Φ𝑖𝑗,2. The classical approach to modeling the PS term, Φ𝑖𝑗, is to split the term (or equivalently the pressure fluctuations term, 𝑝′) into these two distinct contributions of the slow and rapid PS terms, and model separately. Hence, the PS expression of equation (2.30) can be re-written as, (PS) = Φ𝑖𝑗 = Φ𝑖𝑗,1 +Φ𝑖𝑗,2 (2.31) where, Φ𝑖𝑗,1 =14𝜋∫ 𝜌 [𝜕2(𝑢𝑘′ 𝑢𝑙′)𝜕𝑥𝑙𝑥𝑚]∗(𝜕𝑢𝑖′𝜕𝑥𝑗+𝜕𝑢𝑗′𝜕𝑥𝑖)̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ 𝑑 vol|𝐱 − 𝐲| vol (2.32) and, Φ𝑖𝑗,2 =14𝜋∫ 2𝜌 (𝜕𝑈𝑙𝜕𝑥𝑚)∗(𝜕𝑢𝑚′𝜕𝑥𝑙)∗(𝜕𝑢𝑖′𝜕𝑥𝑗+𝜕𝑢𝑗′𝜕𝑥𝑖)̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ 𝑑 vol|𝐱 − 𝐲| vol (2.33) It is worth mentioning that the terms ‘slow’ and ‘rapid’ do not have any physical implications relating to the pace of the process as neither phenomenon precedes the other. Rather, the indication of ‘rapid’ is merely a reference to the fact that the mean velocity gradients (mean strain rate) contribute to the fluctuating pressure immediately, owing to their explicit presence in the equation (2.33). By redistributing the turbulent energy between the Reynolds normal stresses (i = j), the PS term reduces the anisotropy as well as the Reynolds shear stresses (i ≠ j). In the slow PS process, the interactions of the turbulent eddies with one another, results in the reduction of anisotropy. The rapid PS entails a separate physical process, in which the anisotropy is suppressed due to the 50 turbulent fluctuations interacting with the mean velocity gradients (mean flow strain) field. Launder et al. (1975) proposed a model, also referred to as LRR, in which the slow PS term was assumed to be a linear function of the mean velocity gradients, Φ𝑖𝑗,1 = 𝐶1𝜀𝑘(𝜏𝑖𝑗𝑡 +23𝜌𝑘𝛿𝑖𝑗) = −𝐶1𝜌𝜀𝑘(𝑢𝑖′𝑢𝑗′̅̅ ̅̅ ̅̅ −23𝑘𝛿𝑖𝑗) (2.34) Φ𝑖𝑗,2 = −𝐶2 [(𝑃𝑖𝑗 − 𝐶𝑖𝑗) −23𝛿𝑖𝑗(𝑃 − 𝐶)] (2.35) where 𝐶1=1.8, 𝐶2=0.6, 𝑃 =12𝑃𝑘𝑘, and 𝐶 =12𝐶𝑘𝑘. This model is usually referred to as the linear pressure-stress model because of the linear proportionality between Φ𝑖𝑗,1 and 𝑢𝑖′𝑢𝑗′̅̅ ̅̅ ̅̅ , first proposed by Rotta (1951). The model’s coefficients were algebraically, rather than differentially, dependent on the anisotropy of the Reynolds stresses. As previously mentioned, the overall effect of the PS term is to reduce the anisotropy of Reynolds stresses. However, measurements have shown an increase of the anisotropy in the proximity of a solid wall, which has been attributed to the dampening of the turbulent fluctuations in the normal direction to the wall (Versteeg and Malalasekera, 2007). Thus, Launder et al. (1975) introduced additional corrections to equation (2.31), to account for the near-wall effect on the PS model. The additional near wall effect term takes the following form, Φ𝑖𝑗,𝑤 = 𝐶1′𝜀𝑘(𝑢𝑘′ 𝑢𝑚′̅̅ ̅̅ ̅̅ ̅𝑛𝑘𝑛𝑚𝛿𝑖𝑗 −32𝑢𝑖′𝑢𝑘′̅̅ ̅̅ ̅̅ 𝑛𝑗𝑛𝑘 −32𝑢𝑗′𝑢𝑘′̅̅ ̅̅ ̅̅ 𝑛𝑖𝑛𝑘)𝐶𝑙𝑘32⁄𝜀𝑑+ 𝐶2′ (Φ𝑘𝑚,2𝑛𝑘𝑛𝑚𝛿𝑖𝑗 −32Φ𝑖𝑘,2𝑛𝑗𝑛𝑘 −32Φ𝑗𝑘,2𝑛𝑖𝑛𝑘)𝐶𝑙𝑘32⁄𝜀𝑑 (2.36) where 𝐶1′=0.5, 𝐶2′=0.3, 𝑛𝑘 is the component, in the 𝑥𝑘 direction, of the unit normal to the wall, 𝐶𝑙 = 𝐶𝜇3/4 / 𝜅, with 𝜅 being the von Kármán constant (𝜅=0.4187). The LPS model used in this study is based on the correlations proposed by Launder et al. (1975), Gibson and Launder (1978), Fu et al. (1987), Launder (1989), and Launder (1989). The final form for the modelled PS is, 51 (PS) = Φ𝑖𝑗 = Φ𝑖𝑗,1 +Φ𝑖𝑗,2 +Φ𝑖𝑗,𝑤 (2.37) with Φ𝑖𝑗,1, Φ𝑖𝑗,2, and Φ𝑖𝑗,𝑤 defined by equations (2.34) - (2.36). 2.2.2.2.2 Quadratic Pressure-Strain Despite its comprehensive considerations, the LPS model can lead to inaccurate predictions for swirling flows with curved streamlines due to the simple, linear assumption. In order to overcome this deficiency, Speziale et al. (1991) proposed a model that incorporates a quadratic, non-linear correlation for the slow term. This quadratic model, also known as the SSG model, is demonstrated to produce superior performance through more accurate results for highly swirling flows with streamline curvature such as that of hydrocyclones (Cullivan et al., 2003a; Ansys, 2011; Launder and Shima, 1989; Cullivan et al., 2004; Murphy et al., 2007). The quadratic model is implemented in ANSYS Fluent in the following form (Ansys, 2011), Φ𝑖𝑗 = −(𝐶1𝜌𝜀 + 𝐶1∗𝑃)𝑏𝑖𝑗 + 𝐶2𝜌𝜀 (𝑏𝑖𝑘𝑏𝑘𝑗 −13𝑏𝑚𝑛𝑏𝑚𝑛𝛿𝑖𝑗)+ (𝐶3 − 𝐶3∗√𝑏𝑖𝑗𝑏𝑖𝑗) 𝜌𝑘𝑆𝑖𝑗 + 𝐶4𝜌𝑘 (𝑏𝑖𝑘𝑆𝑗𝑘 + 𝑏𝑗𝑘𝑆𝑖𝑘 −23𝑏𝑚𝑛𝑆𝑚𝑛𝛿𝑖𝑗)+ 𝐶5𝜌𝑘(𝑏𝑖𝑘Ω𝑗𝑘 + 𝑏𝑗𝑘Ω𝑖𝑘) (2.38) where 𝑏𝑖𝑗 is the Reynolds anisotropy, and Ω𝑖𝑗 is the mean rotation-rate tensors, which are defined as, 𝑏𝑖𝑗 = −(−𝜌𝑢𝑖′𝑢𝑗′̅̅ ̅̅ ̅̅ +23 𝜌𝑘𝛿𝑖𝑗2𝜌𝑘) (2.39) Ω𝑖𝑗 =12(𝜕𝑈𝑖𝜕𝑥𝑗−𝜕𝑈𝑗𝜕𝑥𝑖) (2.40) The constants take the following default values, 52 𝐶1 = 3.4, 𝐶1∗ = 1.8, 𝐶2 = 4.2, 𝐶3 = 0.8𝐶3∗ = 1.3, 𝐶4 = 1.25, 𝐶5 = 0.4 } (2.41) Note that k is the turbulence kinetic energy and takes the value of the trace of the Reynolds stress tensor as stated in equation (2.17). A transport equation for the turbulence kinetic energy, identical to that of k-ε model, is solved to compute the boundary conditions for the Reynolds stresses. 2.2.2.3 Dissipation Rate The dissipation rate tensor, 𝜀𝑖𝑗, is modelled using Kolmogorov's (1941) hypothesis of local isotropy, 𝜀𝑖𝑗 =23𝛿𝑖𝑗𝜌𝜀 (2.42) This assumption is popular amongst modelers owing to the fact that the dissipative eddies are of the smallest turbulence scales; the dissipation takes place at the smallest scales. Note that the scalar dissipation rate, 𝜀, is modelled by using the same equation as that in the standard k-ε model, i.e. equation (2.26). The six transport equations for the Reynolds stresses, equation (2.10), are solved along with the transport equation for the scalar dissipation rate to close the set of equations. This results in five and seven additional partial differential equations (PDE’s) in 2D and 3D flows, respectively. These additional PDE’s are responsible for the large computational costs incurred by the RSM. 2.2.3 Large Eddy Simulation (LES) LES can be viewed as a compromise between DNS and RANS; the larger scale motion is resolved while the small scale motion is modelled. It would not be incorrect to describe the LES as a model in which DNS is applied to the larger scales and RANS to the smaller scales. Eddies are a characteristic of turbulent flows and have a wide distribution in size and time scale. The larger eddies, which extract energy from the main flow, are anisotropic and significantly influenced by the geometry and boundary conditions. On the contrary, the smaller eddies are nearly isotropic and 53 exhibit a universal behavior. LES computes the larger eddies directly and approximates the smallest eddies, which are accountable for the turbulent kinetic energy dissipation. Therefore, the LES model resolves a larger fraction of scales compared to the RANS models including RSM. However, LES requires finer grids compared to RANS models, which makes it computationally expensive in terms of computer memory (RAM) and CPU time. The essence of LES can be justified as follows (Ferziger, 1996), The larger eddies contain most of the energy The transport of conserved properties is mostly done by the larger eddies The larger eddies are far more flow and geometry dependent; they vary most from flow to flow The smaller eddies show a more universal behavior, are less critical, and are easier to model The distinct feature of the LES model is its spatial filtering operation used to separate the larger scale eddies from the smaller eddies. Instead of ensemble averaging used in RANS, LES employs a filter function, sometimes called filter kernel, and a certain cutoff width. By applying the spatial filtering operation, eddies with length scales greater than the cutoff width are resolved explicitly in an unsteady flow computation process, while the information related to length scales smaller than the cutoff width, are lost. This loss of small scale information, as well as the interaction between the larger scale and smaller scale eddies, makes it necessary to model the small eddies effect, referred to as subgrid scale (SGS) stresses. The following subsections provide more details about the LES model and SGS stresses. 2.2.3.1 Filtered Navier-Stokes Equations ?̅?(𝒙, 𝑡) = ∫ ∫ ∫ 𝐺(𝒙, 𝒙′, Δ)𝜙(𝒙′, 𝑡) 𝑑𝑥1′ 𝑑𝑥2′ 𝑑𝑥3′+∞−∞+∞−∞+∞−∞ (2.43) where ?̅?(𝒙, 𝑡) ≡ filtered function and 𝜙(𝒙, 𝑡) ≡ original, unfiltered function and Δ ≡ filter cutoff width 54 The most common form of the filtering function, or filter kernel, for finite volume method (FVM) in 3D LES simulations, is the top-hat or volume-average box filter (Deardorff, 1970), defined as, 𝐺(𝒙, 𝒙′, ∆) = 𝑓(𝑥) = {1/∆3 |𝒙 − 𝒙′| ≤ ∆/20 |𝒙 − 𝒙′| > ∆/2 (2.44) Many other filters have been proposed, amongst which the Gaussian filter and the spectral cutoff filter (Fourier cutoff filter) are more commonly used. Refer to the following turbulence texts for a more in-depth discussion (Ferziger, 1996; Sagaut, 2001; Wilcox, 2006; Versteeg and Malalasekera, 2007). The cutoff width, Δ, provides an indicative size criterion for eddies used to determine whether to retain (and resolve) them in the calculations ( > Δ), or to reject (and model) them. Δ can, in principle, be selected to be any size but it is sensible to choose a cutoff width no smaller than the grid size in CFD calculations with the finite volume method (FVM); details finer than the grid size are lost anyways. The most widely employed choice is for the cutoff width to be of the same order as the grid size. Therefore, in three-dimensional computations (including this study) the cutoff width is selected as, ∆= (∆𝑥 ∆𝑦 ∆𝑧)1 3⁄ (2.45) with Δx, Δy, and Δz being grid cells dimensions. Combining equations (2.43) - (2.45), results in the following filtered form for the FVM used in this study, 𝜙?̅?(𝒙, 𝑡) =1𝑉∫ ∫ ∫ 𝜙𝑖(𝝃, 𝑡) 𝑑𝜉 𝑑𝜂 𝑑𝜁𝑥+12∆𝑥𝑥−12∆𝑥𝑦+12∆𝑦𝑦−12∆𝑦𝑧+12∆𝑧𝑧−12∆𝑧 (2.46) where 𝜙?̅? is the resolvable-scale filtered quantity and 𝑉 = Δ3 is the cell volume. It is worth mentioning that the subgrid-scale (SGS) quantity, 𝜙𝑖′, is defined as, 𝜙𝑖′(𝒙, 𝑡) = 𝜙𝑖(𝒙, 𝑡) − ?̅?(𝒙, 𝑡) (2.47) Note the similarity between the SGS quantities in LES and the fluctuating quantities in the RSM (and RANS in general). The filter kernel, applied here, is, 55 𝐺(𝒙, 𝒙′, ∆) = {1∆3=1𝑉 𝒙′ ∈ 𝑉0 𝒙′ ∉ 𝑉 (2.48) where the condition 𝒙′ ∈ 𝑉 means the point of interest, with the coordinates 𝒙′, belongs to the computational cell V. Applying the above-mentioned filter kernel to equations (2.1) and (2.2) leads to the filtered unsteady Navier-Stokes equations, 𝜕𝑢?̅?𝜕𝑥𝑖= 0 (2.49) 𝜕(𝜌𝑢?̅?)𝜕𝑡+𝜕𝜕𝑥𝑗(𝜌𝑢𝑖𝑢𝑗̅̅ ̅̅ ̅) = −𝜕?̅?𝜕𝑥𝑖+ 𝜇𝜕2𝑢?̅?𝜕𝑥𝑗𝜕𝑥𝑗 (2.50) where the quantities with the bar are filtered (not time-averaged) quantities. All the terms in equation (2.50) are exact except for the convective flux, i.e. 𝑢𝑖𝑢𝑗̅̅ ̅̅ ̅. Therefore, 𝜌𝑢𝑖𝑢𝑗̅̅ ̅̅ ̅ ≠ 𝜌𝑢?̅?𝑢?̅? (2.51) A correlation must be introduced to approximate the difference between the two sides of the inequality in equation (2.51), which is called subgrid-scale (SGS) stresses, 𝜏𝑖𝑗, 𝜏𝑖𝑗 = 𝜌𝑢𝑖𝑢𝑗̅̅ ̅̅ ̅ − 𝜌𝑢?̅?𝑢?̅? (2.52) thus rendering the filtered equation (2.50) solvable computationally. Using the decomposition of equation (2.47) on the convective flux in equation (2.50) leads to, 𝜌𝑢𝑖𝑢𝑗̅̅ ̅̅ ̅ = 𝜌(𝑢?̅? + 𝑢𝑖′)(𝑢?̅? + 𝑢𝑗′)̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅ = 𝜌𝑢?̅?𝑢?̅?̅̅ ̅̅ ̅ + 𝜌𝑢?̅?𝑢𝑗′̅̅ ̅̅ ̅ + 𝜌𝑢𝑖′𝑢?̅?̅̅ ̅̅ ̅ + 𝜌𝑢𝑖′𝑢𝑗′̅̅ ̅̅ ̅̅ = 𝜌𝑢?̅?𝑢?̅? + (𝜌𝑢?̅?𝑢?̅?̅̅ ̅̅ ̅ − 𝜌𝑢?̅?𝑢?̅?) + 𝜌𝑢?̅?𝑢𝑗′̅̅ ̅̅ ̅ + 𝜌𝑢𝑖′𝑢?̅?̅̅ ̅̅ ̅ + 𝜌𝑢𝑖′𝑢𝑗′̅̅ ̅̅ ̅̅ (2.53) Substituting equation (2.53) into equation (2.52), the SGS stresses can be expressed as, 56 𝜏𝑖𝑗 = (𝜌𝑢?̅?𝑢?̅?̅̅ ̅̅ ̅ − 𝜌𝑢?̅?𝑢?̅?)⏟ 𝐿𝑖𝑗+ (𝜌𝑢?̅?𝑢𝑗′̅̅ ̅̅ ̅ + 𝜌𝑢𝑖′𝑢?̅?̅̅ ̅̅ ̅)⏟ 𝐶𝑖𝑗+ 𝜌𝑢𝑖′𝑢𝑗′̅̅ ̅̅ ̅̅⏟ 𝑅𝑖𝑗 (2.54) where 𝐿𝑖𝑗, 𝐶𝑖𝑗, 𝑅𝑖𝑗 are known as Leonard stresses, cross term (stresses), and SGS (or LES) Reynolds stresses. The physical significance of these terms is described as follows, Leonard stresses, also called outscatter term, represent the interaction between resolved eddies to produce small scale eddies, and can be calculated explicitly from the filtered velocity field and do not require modelling (Leonard, 1975). This term exists due to the fact that in space filtering (except for the Fourier cutoff filter), unlike time or ensemble averaging, a second filtering operation changes the filtered quantity, i.e. ?̅̅? ≠ ?̅?. Leonard stresses have been shown to remove a significant amount of energy from the resolvable scale. Cross stresses (cross term) is due to the interaction between resolved scale eddies and small (SGS) eddies. These interactions can result in energy transfer in either direction, although on average this term removes energy from the resolved scale to the small scale. LES Reynolds stresses, also known as SGS Reynolds stresses and true subgrid scale term, represent convective momentum transfer due the interaction between small (SGS) eddies. Since it transfers the produced energy from the small scale interactions to the large, resolvable scale, it is also called the backscatter term. Equation (2.50) can now be re-arranged as, 𝜕(𝜌𝑢?̅?)𝜕𝑡+𝜕𝜕𝑥𝑗(𝜌𝑢?̅?𝑢?̅?) = −𝜕?̅?𝜕𝑥𝑖+ 𝜇𝜕2𝑢?̅?𝜕𝑥𝑗𝜕𝑥𝑗−𝜕𝜏𝑖𝑗𝜕𝑥𝑗 (2.55) Note the similarity between the LES momentum equation (2.55) and the RANS equation (2.9); the difference is that the overhead bar in LES indicates filtered variables, whereas in RANS it represents time-averaged variables. Moreover, 𝜏𝑖𝑗 in LES represents SGS stresses, whereas it represents Reynolds stresses in RANS. The role SGS stresses play in the LES momentum equation is similar to the role Reynolds stresses play in RANS. As was the case with RANS, the fundamental problem of LES is to find a suitable approximation for 𝜏𝑖𝑗, to accurately model SGS stresses. This 57 study employs the dynamic procedure (model) which makes the isotropic turbulence assumption for the small scale eddies. The basis of the dynamic model is the Smagorinsky-Lilly model. The following subsections present more detail about both the Smagorinsky-Lilly model and dynamic procedure. 2.2.3.2 Smagorinsky-Lilly SGS Model The most widely used subgrid scale model was proposed by Smagorinsky (1963) in which the Boussinesq hypothesis was adopted with the following justification; since the smallest scale eddies are almost isotropic, anisotropy declines considerably with the length scale, this hypothesis should provide a reasonably accurate description of the unresolved eddies, i.e. SGS stresses. Therefore, the Smagorinsky model is essentially an eddy viscosity model which assumes that the SGS stresses are proportional to the strain rate for the resolved flow (similar to the proportionality between the Reynolds stresses and the strain rate in RANS), 𝜏𝑖𝑗 −13𝜏𝑘𝑘𝛿𝑖𝑗 = −2𝜇𝑡 𝑆 𝑖𝑗 (2.56) where the overbar indicates spatial averaging. Thus, the resolved strain tensor is, 𝑆 𝑖𝑗 = 12(𝜕𝑢𝑖𝜕𝑥𝑗+𝜕𝑢𝑗𝜕𝑥𝑖) (2.57) Note that the second term on the LHS ensures that the sum of normal SGS stresses is equal to the SGC kinetic energy. Similar to the turbulent viscosity defined for the RANS equations, the SGS turbulent viscosity can be described in terms of a length scale and a velocity scale. The filter cutoff width, Δ, (or alternatively Cs Δ), would be a straightforward choice for the length scale as it essentially determines the SGS eddies. The velocity scale, similar to RANS, is stated as the product of the length scale and the (spatial) average of the (resolved) strain rate. For example, ∆ × |𝑆̅| (or 𝐶𝑠∆ × |𝑆̅| ), with |𝑆̅| = √2𝑆̅ 𝑖𝑗 𝑆̅ 𝑖𝑗 . The subgrid length scale, 𝐿𝑠 , employed in this study is determined from the minimum of the following two values, 58 𝐿𝑠 = 𝑚𝑖𝑛(𝜅𝑑, 𝐶𝑠∆) (2.58) where 𝜅, 𝑑, 𝐶𝑠, 𝑎𝑛𝑑 Δ are the von Kármán constant, the distance to the closest wall, the Smagorinsky constant, or parameter, and the filter cutoff width, respectively. The filter cutoff width (or local grid scale), Δ, is simply the cube root of the volume of the computational cell, i.e. Δ = 𝑉1 3⁄ . Hence, the subgrid-scale turbulent eddy viscosity, 𝜇𝑡 (or 𝜇𝑆𝐺𝑆 as denoted in some texts), is computed as, 𝜇𝑡 = 𝜌𝐿𝑠2 |𝑆| (2.59) It is worth mentioning that Smagorinsky constant 𝐶𝑠 is not a universal constant. Nevertheless, a value of ≈ 0.1 has been shown to produce satisfactory results for a broad range of flows (Smagorinsky 1963). 2.2.3.3 Dynamic Model Germano et al. (1991) proposed a model, later improved by Lilly (1992), to compute the Smagorinsky parameter in a dynamic way, by using the information from the resolved velocity field as a priori estimation of the SGS model parameter. Therefore, the model constant need not be specified in advance, rather it is computed at every spatial point and time step dynamically. This is why the model is referred to as dynamic model, even though strictly speaking it is a procedure and not a model since its basis is a previously developed model, the Smagorinsky model. The concept of the dynamic procedure is that a second filter, usually called test filter, is applied to the Navier-Stokes equations. The test filter width, ∆̃, is twice as large as the grid filter width, Δ. The two filtering procedures produce their respective resolved flow fields. The difference between these resolved fields represents the contribution of the small scales with their size lying between the grid filter and the test filter, i.e. ∆ < 𝑠𝑚𝑎𝑙𝑙 𝑠𝑐𝑎𝑙𝑒𝑠 𝑠𝑖𝑧𝑒 < ∆̃. The information obtained from these scales is what the dynamic model uses to evaluate the model’s constant. The following presents a mathematical description of the dynamic procedure. The subgrid scale stresses in the test-filter LES can be written as (~ sign indicates test filter variables), 59 𝑇𝑖𝑗 = 𝜌𝑢𝑖𝑢𝑗̅̅ ̅̅ ̅̃ − 𝜌𝑢?̃̅?𝑢?̃̅? (2.60) The large scale contributions of the SGS stresses at the test filter level, which can be computed directly, is in fact the Leonard stresses related to the test-filter level, 𝐿𝑖𝑗 = 𝜌𝑢?̅?𝑢?̅?̃ − 𝜌𝑢?̃̅?𝑢?̃̅? (2.61) It can be seen from equations (2.52), (2.61), and (2.62) that, 𝐿𝑖𝑗 = 𝑇𝑖𝑗 − ?̃? 𝑖𝑗 (2.62) The above equation describes a mathematical identity, called Germano’s identity, which gives the basis of the dynamic model. When the Smagorinsky model is used as a base model for the dynamic procedure, Germano’s identity yields, 𝐿𝑖𝑗 −13𝐿𝑘𝑘𝛿𝑖𝑗 = 𝐶𝑀𝑖𝑗 (2.63) with 𝑀𝑖𝑗 = −2𝜌 (∆̃2|𝑆̅̃|𝑆̅̃ − ∆2|𝑆̅|𝑆̅ 𝑖𝑗 ̃ ) (2.64) The constant C is computed based on the least square method suggested by Lilly (1992), 𝐶 = 𝐶𝑠2 =〈𝐿𝑖𝑗𝑀𝑖𝑗〉〈𝑀𝑖𝑗𝑀𝑖𝑗〉 (2.65) where 〈 〉 indicates a spatial averaging over the domain. More details on these models and the numerical calculations can be found in turbulence books by Ferziger (1996) and Sagaut (2001). 2.3 Three-Dimensional Modeling The RANS equations for unsteady, three-dimensional flows with constant properties can be expressed in the cylindrical coordinates, in the vectorial form, as follows, 60 𝜌𝜕𝑈𝜕𝑡+ 𝜌(?⃗? . ∇⃗ )?⃗? = −∇⃗ 𝑃 + 𝜇∇2?⃗? − ∇⃗ . 𝜏?̿? (2.66) where the velocity vector in the cylindrical coordinate is ?⃗? (𝑟, 𝜃, 𝑧) = 𝑢𝑟(𝑟, 𝜃, 𝑧)?̂?𝑟 +𝑢𝜃(𝑟, 𝜃, 𝑧)?̂?𝜃 + 𝑢𝑧(𝑟, 𝜃, 𝑧)?̂?𝑧 with ?̂?𝑟, ?̂?𝜃, and ?̂?𝑧 the unit normal vectors in their respective directions of radial, tangential, and axial. The Reynolds stresses tensor, in the cylindrical coordinate, takes the following form, 𝜏?̿? = [𝜏𝑟𝑟𝑡 𝜏𝑟𝜃𝑡 𝜏𝑟𝑧𝑡𝜏𝜃𝑟𝑡 𝜏𝜃𝜃𝑡 𝜏𝜃𝑧𝑡𝜏𝑧𝑟𝑡 𝜏𝑧𝜃𝑡 𝜏𝑧𝑧𝑡] = −𝜌 [𝑢𝑟′𝑢𝑟′̅̅ ̅̅ ̅̅ 𝑢𝑟′𝑢𝜃′̅̅ ̅̅ ̅̅ ̅ 𝑢𝑟′𝑢𝑧′̅̅ ̅̅ ̅̅𝑢𝜃′ 𝑢𝑟′̅̅ ̅̅ ̅̅ ̅ 𝑢𝜃′ 𝑢𝜃′̅̅ ̅̅ ̅̅ ̅ 𝑢𝜃′ 𝑢𝑧′̅̅ ̅̅ ̅̅𝑢𝑧′𝑢𝑟′̅̅ ̅̅ ̅̅ 𝑢𝑧′𝑢𝜃′̅̅ ̅̅ ̅̅ 𝑢𝑧′𝑢𝑧′̅̅ ̅̅ ̅̅] (2.67) It is worth mentioning that the Reynolds stresses can be converted to the Cartesian coordinate notation, using the following conversion matrix, 𝜏?̿? = [𝜏𝑟𝑟𝑡 𝜏𝑟𝜃𝑡 𝜏𝑟𝑧𝑡𝜏𝜃𝑟𝑡 𝜏𝜃𝜃𝑡 𝜏𝜃𝑧𝑡𝜏𝑧𝑟𝑡 𝜏𝑧𝜃𝑡 𝜏𝑧𝑧𝑡]= [𝑐𝑜𝑠𝜃 𝑠𝑖𝑛𝜃 0−𝑠𝑖𝑛𝜃 𝑐𝑜𝑠𝜃 00 0 1] [𝜏𝑥𝑥𝑡 𝜏𝑥𝑦𝑡 𝜏𝑥𝑧𝑡𝜏𝑦𝑥𝑡 𝜏𝑦𝑦𝑡 𝜏𝑦𝑧𝑡𝜏𝑧𝑥𝑡 𝜏𝑧𝑥𝑡 𝜏𝑧𝑧𝑡] [𝑐𝑜𝑠𝜃 −𝑠𝑖𝑛𝜃 0𝑠𝑖𝑛𝜃 𝑐𝑜𝑠𝜃 00 0 1] The terms in equation (2.66) are expanded to show their three-dimensional form in the cylindrical coordinate. The convective derivative of the LHS has the following expanded form in terms of the velocity components, (2.68) (?⃗? . ∇⃗ )?⃗? =[ 𝑈𝑟𝜕𝑈𝑟𝜕𝑟+𝑈𝜃𝑟𝜕𝑈𝑟𝜕𝜃−𝑈𝜃𝑈𝜃𝑟+ 𝑈𝑧𝜕𝑈𝑟𝜕𝑧𝑈𝑟𝜕𝑈𝜃𝜕𝑟+𝑈𝜃𝑟𝜕𝑈𝜃𝜕𝜃+𝑈𝑟𝑈𝜃𝑟+ 𝑈𝑧𝜕𝑈𝜃𝜕𝑧𝑈𝑟𝜕𝑈𝑧𝜕𝑟+𝑈𝜃𝑟𝜕𝑈𝑧𝜕𝜃+ 𝑈𝑧𝜕𝑈𝑧𝜕𝑧 ] (2.69) The pressure gradient is, 61 ∇⃗ 𝑃 =[ 𝜕𝑃𝜕𝑟1𝑟𝜕𝑃𝜕𝜃𝜕𝑃𝜕𝑧 ] (2.70) The Laplacian in the viscous forces term takes the form as below, ∇2?⃗? =[ 𝜕2𝑈𝑟𝜕𝑟2+1𝑟2𝜕2𝑈𝑟𝜕𝜃2+𝜕2𝑈𝑟𝜕𝑧2+1𝑟𝜕𝑈𝑟𝜕𝑟−2𝑟2𝜕𝑈𝜃𝜕𝜃−𝑈𝑟𝑟2𝜕2𝑈𝜃𝜕𝑟2+1𝑟𝜕𝑈𝜃𝜕𝑟+1𝑟2𝜕2𝑈𝜃𝜕𝜃2+2𝑟2𝜕𝑈𝑟𝜕𝜃−𝑈𝜃𝑟2+𝜕2𝑈𝜃𝜕𝑧2𝜕2𝑈𝑧𝜕𝑟2+1𝑟𝜕𝑈𝑧𝜕𝑟+1𝑟2𝜕2𝑈𝑧𝜕𝜃2+𝜕2𝑈𝑧𝜕𝑧2 ] (2.71) Finally, the divergence of the Reynolds stresses assumes the following form, ∇⃗ . 𝜏̿ =[ 𝜕𝜏𝑟𝑟𝜕𝑟+1𝑟(𝜕𝜏𝜃𝑟𝜕𝜃+ 𝜎𝑟𝑟 − 𝜎𝜃𝜃) +𝜕𝜏𝑧𝑟𝜕𝑧𝜕𝜏𝑟𝜃𝜕𝑟+1𝑟(𝜕𝜏𝜃𝜃𝜕𝜃+ 𝜏𝑟𝜃 + 𝜏𝜃𝑟) +𝜕𝜏𝑧𝜃𝜕𝑧𝜕𝜏𝑟𝑧𝜕𝑟+1𝑟(𝜕𝜏𝜃𝑧𝜕𝜃+ 𝜏𝑟𝑧) +𝜕𝜏𝑧𝑧𝜕𝑧 ] (2.72) 2.4 Axisymmetric Modeling The RANS equations in the axisymmetric model are taken from their three-dimensional counterparts, with an additional assumption that changes in the tangential direction are negligible; therefore, rendering derivatives in the tangential direction equal to zero, i.e. 𝜕 𝜕𝜃⁄ =0. Hence, the RANS equations in the axisymmetric model are, Continuity: 𝜕𝑈𝑧𝜕𝑧+𝜕𝑈𝑟𝜕𝑟+𝑈𝑟𝑟= 0 (2.73) 62 Radial momentum: 𝜌 (𝜕𝑈𝑟𝜕𝑡+ 𝑈𝑟𝜕𝑈𝑟𝜕𝑟−𝑈𝜃2𝑟+ 𝑈𝑧𝜕𝑈𝑟𝜕𝑧)= −𝜕𝑃𝜕𝑟+ 𝜇 [1𝑟𝜕𝜕𝑟(𝑟𝜕𝑈𝑟𝜕𝑟) −𝑈𝑟𝑟2+𝜕2𝑈𝑟𝜕𝑧2] − 𝜌𝜕𝜕𝑟(𝑢′𝑟2)−𝜌𝑢′𝑟2 − 𝜌𝑢′𝜃2𝑟−𝜕𝜕𝑧(𝜌𝑢′𝑧𝑢′𝑟) (2.74) Tangential momentum: 𝜌 (𝜕𝑈𝜃𝜕𝑡+ 𝑈𝑟𝜕𝑈𝜃𝜕𝑟+𝑈𝑟𝑈𝜃𝑟+ 𝑈𝑧𝜕𝑈𝜃𝜕𝑧)= 𝜇 [1𝑟(𝜕𝜕𝑟(𝑟𝜕𝑈𝜃𝜕𝑟)) −𝑈𝜃𝑟2+𝜕2𝑈𝜃𝜕𝑧2] −𝜕𝜕𝑟(𝜌𝑢′𝑟𝑢′𝜃) − 2𝜌𝑢′𝑟𝑢′𝜃𝑟−𝜕𝜕𝑧(𝜌𝑢′𝑧𝑢′𝜃) (2.75) Axial momentum: 𝜌 (𝜕𝑈𝑧𝜕𝑡+ 𝑈𝑟𝜕𝑈𝑧𝜕𝑟+ 𝑈𝑧𝜕𝑈𝑧𝜕𝑧)= −𝜕𝑃𝜕𝑧+ 𝜇 [1𝑟𝜕𝜕𝑟(𝑟𝜕𝑈𝑧𝜕𝑟) +𝜕2𝑈𝑧𝜕𝑧2] −𝜕𝜕𝑟(𝜌𝑢𝑧′𝑢𝑟′̅̅ ̅̅ ̅̅ ) −𝜌𝑢𝑧′𝑢𝑟′̅̅ ̅̅ ̅̅ ̅̅𝑟−𝜕𝜕𝑧(𝜌𝑢′𝑧2) (2.76) 2.5 Ideal Vortex Laws The two ideal cases for swirling flows are the free-vortex and forced-vortex (solid-body rotation) flows. The tangential velocity profiles for these two ideal vortex flows are derived in the following subsections from the θ-direction momentum equation; equation (2.2) for θ direction in the cylindrical coordinate is written as, 63 𝜌 (𝜕𝑢𝜃𝜕𝑡+ 𝑢𝑟𝜕𝑢𝜃𝜕𝑟+𝑢𝜃𝑟𝜕𝑢𝜃𝜕𝜃+𝑢𝑟𝑢𝜃𝑟+ 𝑢𝑧𝜕𝑢𝜃𝜕𝑧)= −1𝑟𝜕𝑝𝜕𝜃− [𝜕𝜏𝑟𝜃𝜕𝑟+1𝑟(𝜕𝜏𝜃𝜃𝜕𝜃+ 𝜏𝑟𝜃 + 𝜏𝜃𝑟) +𝜕𝜏𝑧𝜃𝜕𝑧] (2.77) 2.5.1 Forced Vortex Equation (2.77) can be greatly simplified for steady, incompressible, axisymmetric vortex flow (𝜕 𝜕𝜃⁄ =0), with no velocity in the radial and axial directions, i.e. 𝑢𝑟 = 𝑢𝑧 = 0, to, 𝜕𝜏𝑟𝜃𝜕𝑟+ 2𝜏𝑟𝜃𝑟=1𝑟2𝜕𝜕𝑟(𝑟2𝜏𝑟𝜃) = 0 (2.78) The viscous stress component (expressed in equation (2.5) for the Cartesian coordinates) takes the following form in the cylindrical coordinate, using the constant Newtonian viscosity hypothesis, 𝜏𝑟𝜃 = 𝜇 [𝑟𝜕𝜕𝑟(𝑢𝜃𝑟) +1𝑟𝜕𝑢𝑟𝜕𝜃] (2.79) Note that the second term in the viscous stress component vanishes as there is no gradient in the θ direction. Substituting equation (2.79) into equation (2.78) leads to, 𝜕𝜕𝑟[𝜇𝑟3𝜕𝜕𝑟(𝑢𝜃𝑟)] = 0 𝜇≠0 ⇒ 𝜕𝜕𝑟[𝑟3𝜕𝜕𝑟(𝑢𝜃𝑟)] = 0 (2.80) A trivial solution is when μ=0; any velocity profile will satisfy the above equation. Since the derivative is equal to zero, the term inside the bracket must be equal to a constant, 𝑟3𝜕𝜕𝑟(𝑢𝜃𝑟) = −𝐶0 ⟹ 𝑢𝜃 =𝐶0𝑟+ 𝐶1𝑟 (2.81) No singularity condition (𝑢𝜃 ≠ 0 @ 𝑟 = 0) dictates that 𝐶0 = 0, therefore, 𝑢𝜃 = 𝐶1𝑟 = 𝜔𝑟 (2.82) 64 where the constant ω is angular velocity. The above equation states a solid-body-rotation velocity profile. This type of motion is an ideal vortex flow and is referred to as forced-vortex flow, solid-body-rotation flow, or rigid-body-rotation flow. 2.5.2 Free Vortex The other ideal vortex motion is when the fluid flow is inviscid, i.e. 𝜏𝑟𝜃 = 𝜏𝑧𝜃 = 𝜏𝑟𝑧 = 0. Thus, equation (2.77) is simplified to, 𝑢𝑟 (𝜕𝑢𝜃𝜕𝑟+𝑢𝜃𝑟) = 0 (2.83) Note that to arrive at this equation, radial motion of fluid elements were allowed (𝑢𝑟 ≠ 0). The solution to equation (2.83) is achieved by the separation of variables, 𝑢𝜃 =𝐾𝑟 (2.84) where K is a constant. This equation states the tangential velocity profile for the free vortex flow (also called loss-free vortex). The moment of momentum of fluid elements in the radial direction is conserved in this ideal type of vortex flow. As discussed in Chapter 1, the tangential velocity profile in the hydrocyclone is a combination of these two types of the ideal vortex flows.65 Chapter 3: Effect of Swirl on Turbulence This chapter presents CFD simulations of the turbulent swirling flow inside a rotating pipe (D=2cm and L/D=20), conducted using the quadratic pressure-strain Reynolds Stress Model. The results give insight into the complex hydrodynamics of turbulent swirling flows with potential applications in swirling separation equipment including the hydrocyclone. A by-product of this understanding will be the introduction of a novel criterion which would enable researchers and engineers to determine whether complex and expensive turbulence models, such as RSM, could be avoided, without compromising accuracy, in favour of simpler and less computationally-intensive models such as the k-ε model. In order to better understand the interaction between forces, and to simplify the cyclone geometry, the response of the turbulent pipe flow to rotation of the pipe wall was examined. Insights into this type of problem can be “very useful for improving turbulence model related to swirling flows” and enhance our understanding of the physics involved (Orlandi and Fatica, 1997). The present study uses this model problem, which is well established in the literature, to ascertain the conditions under which computationally cheaper turbulence models such as k-ε model should be accurate. These conditions will provide guidelines for the choice of turbulence model required for desired accuracy in modeling the flow field within more complex geometries such as the hydrocyclone. This will enable us to make a proper and economical choice of turbulence model as a design optimization tool. The simulations were carried out for four different Reynolds numbers and a range of rotation numbers (N). The present study adds a further contribution to the literature by introducing a novel integral measure for analyzing the effect of rotation on turbulence. The main objective of the study is to shed light on the response of and the interplay between the convective, pressure, and turbulence forces to changes in rotation number and Reynolds number. Additionally, the competition between the forces in establishing the flow throughout the domain is evaluated. It should be noted that due to limited computational facilities and the relatively high computational costs incurred by RSM, the pipe length used in the simulations is primarily restricted to L = 20 D. 66 3.1 Numerical Scheme and Boundary Conditions Commercial CFD package ANSYS Fluent was employed to perform the calculations for the unsteady, axisymmetric flow of water inside a pipe with the dimeter of D=2 cm and length of L= 40 cm (L/D=20). The simulations were carried out for four different axial Reynolds numbers of 3500, 7000, 14000 and 28000 and eight rotation numbers, N = 0 (stationary pipe), 0.14, 0.29, 0.57, 0.87, 1.14, 1.43, and 2.86, making a total of 32 simulation cases. The Semi-Implicit Method for Pressure Linked Equations (SIMPLE) was used for pressure-velocity coupling. The Pressure Staggering Option (PRESTO!), recommended for highly swirling flows (Ghadirian et al., 2013), was used for pressure spatial discretization. The QUICK algorithm (Quadratic Upwind Interpolation for Convective Kinematics) was employed for spatial discretization of the momentum, swirl and Reynolds stresses equations. This high order scheme was employed to minimize numerical diffusion and improve convergence (Ghadirian et al., 2013). The inlet boundary condition was assigned velocity inlet whereas the outlet boundary condition was set to pressure outlet. The inlet was assigned uniform velocity condition whereas pressure outlet was chosen for the outlet condition. The hydraulic diameter was defined as the pipe diameter with a turbulent intensity of 10%. The k-ε model was used for initialization, before switching to the linear pressure-strain RSM to march in time for additional time steps. Once a crude flow field was established, i.e. pressure and momentum were in balance, the turbulence model was switched to the quadratic pressure-strain RSM. This initialization process was used to avoid numerical instability and convergence difficulties caused by the QPS model. An overall time step size of 10-4 to 10-3 s was used. The average computational time for each case was approximately between 240 - 336 hours, using 12 cores on an i7 CPU machine. The number of mesh elements was 960,000. The evolution of several flow variables, including the velocity components and their derivatives, wall shear stress, and the pressure gradient, were monitored for the convergence assessment. Simulations were deemed converged when the monitored variables leveled off. The converged solutions exhibited no considerable changes (less than 1% in most cases), neither with increasing the time step size nor after numerous time steps (minimum 1000). 67 3.2 Axisymmetric vs. Three-Dimensional Simulations In order to verify the choice of axisymmetric model, a three-dimensional (3D) simulation was carried out on a fine grid with 4.6 million cells and identical boundary and flow conditions. A comparison between the radial profiles of non-dimensionalized axial and swirl velocity have been presented in Figure 3.1 and Figure 3.2. It can be seen that the velocity field predictions of the axisymmetric and 3D simulations nearly perfectly coincide. Therefore, the simulations were carried out using the axisymmetric model to reduce the computational cost, without compromising accuracy. Figure 3.1: Dimensionless axial velocity at x/D=20 at Re=7000 and N=0.29; normalized by the bulk axial velocity. 68 Figure 3.2: Dimensionless swirl velocity at x/D=20 at Re=7000 and N=0.29; normalized by the pipe wall circumferential velocity. 3.3 Grid Sensitivity The grid sensitivity study was performed using three different grids, namely coarse, intermediate and fine, with 60,000, 240,000, and 960,000 mesh elements, respectively. Figure 3.3 compares between the axial velocity obtained from these grids for the condition of Re = 7000 and N = 0. It can be seen that the predictions of the intermediate and fine grids collapse onto one another. A similar trend was observed for the rotating pipe at N = 1.14. The simulations were carried out using the fine grid to ensure maximum accuracy. 69 Figure 3.3: Axial velocity profiles for three different grids at x / D =20 at Re = 7000 and N = 0. 3.4 Experimental Validation Figure 3.4 compares the axial velocity profiles obtained from the present work with the RSM simulations of Hirai et al. (1988) and the experimental data of Kikuyama et al. (1983) for two different rotation numbers of N = 0.5 and N = 1. It can be seen that the results of the present simulations are in reasonably good agreement with both sets of data. 70 Figure 3.4: Comparison between the normalized axial velocity profiles from the present work (PW), Hirai et al.’s RSM results, and Kikuyama et al.’s experimental data at Re = 10000 for two rotation numbers N = 0.5 and N = 1. 3.5 Pressure Field The pressure field evolution, along the axial direction, is presented in this section. The pressure difference between planes of interest and the inlet were normalized by the dynamic pressure of the mean flow, i.e. 𝜌𝑈𝑚2 2⁄ , forming the Euler number as follows, 𝐸𝑢 =𝑝 − 𝑝𝑖𝑛12𝜌𝑈𝑚2 (3.1) where p is the area-weighted average of pressure at the plane of interest and pin is the area-weighted average of pressure at the inlet. The dimensionless form of the pressure profile (Eu) has been illustrated in Figure 3.5 for Re = 7000 and a range of pipe wall rotational speeds. It was observed that the Euler number, i.e. total pressure loss along the pipe, decreases with increasing rotational speed which is attributed to the laminarization phenomenon (Murakami and Kikuyama, 1980). Note that the pressure drop in the immediate downstream of the inlet (x/D ≤ 5) remains rather unaffected by increasing the rotation numbers to 1.43, corresponding to a rotational speed of 0.5 rad/s. 71 Figure 3.5: Euler number distribution along the dimensionless pipe length, normalized by diameter, at Re=7000 and a range of rotation numbers. 3.6 Velocity Field Non-dimensionalized axial velocity profiles at different axial positions, for Re = 7000 and N = 1.43, are illustrated in Figure 3.6. The laminarization of the flow field is evident from the evolution of the boundary layer, from turbulent at planes near the inlet into a laminar-like profile at distances farther away from the inlet. Moreover, it can be seen that the velocity profiles at x = 18D and x = 19D collapse onto each other, indicating that the flow has reached its fully developed state. Figure 3.7 shows the non-dimensionalized swirl velocity at different x/D for Re = 7000 and N = 1.43. It can be seen that the swirl velocity field is almost fully established at x >18D, with insignificant changes between the swirl velocity beyond x = 18D. These results are similar to the observation for the axial velocity field presented in Figure 3.6. This represents the time required for the swirling turbulent boundary layer to establish itself for the given flow conditions. Moreover, it is interesting to note that at x/D = 4 the influence of the rotational wall is mostly limited to a radial distance of 0.4R, i.e. 0.6 < r/R < 1, from the rotating wall. Further downstream at x/D = 18 this influence extends to almost 70% of the radial domain away from the rotating wall (0.3 ≤ r/R ≤ 1). This suggests that the central core of the flow (r/R ≤ 0.3) will not be affected 72 significantly with an increase in the pipe length. However, further increase of the pipe rotational speed would be expected to extend the reach of the swirl influence to the flow central core. Figure 3.6: Axial velocity profiles, normalized by mean axial velocity, at different axial locations for Re=7 000 and N=1.43. Figure 3.7: Swirl velocity profiles, normalized by the pipe wall rotational velocity, at different axial locations for Re=7 000 and N = 1.43. 73 The laminarization phenomenon is attributed to the suppression of turbulence, e.g. −𝜌𝑢𝑧′𝑢𝑟′̅̅ ̅̅ ̅̅ , due to the swirl (Hira et al., 1988). Figure 3.8 illustrates the profiles of the Reynolds shear stress,−𝜌𝑢𝑧′𝑢𝑟′̅̅ ̅̅ ̅̅ , at Re = 7000 for a variety of rotation numbers. It is evident that an increase in the rotation number is followed by a decrease in the Reynolds stress, −𝜌𝑢𝑧′𝑢𝑟′̅̅ ̅̅ ̅̅ . In other words, increasing the swirl leads to dampening the turbulence and therefore diminishes the turbulent momentum exchange. Note that at higher rotation numbers, the Reynold shear stress profile flattens out, indicating the turbulence suppression. Note that the Reynolds shear stress, at all N’s, reaches its maximum in the close proximity of the wall, due to the strong convective momentum transfer between the pipe wall and the flow in its vicinity. Figure 3.8: Reynolds shear stress at Re = 7 000 and x = 19D for different rotation numbers. 3.7 Force Balance Analysis A force balance analysis has been performed by numerically calculating the different terms of the radial momentum balance, i.e. equation (2.74), to quantify the contribution of each term to the total force balance. The sum of the terms on the LHS of equation (2.74) forms the convective forces, whereas the pressure gradient term is denoted by pressure force. Viscous forces are represented by the second term on the right-hand side (RHS), whereas the summation of the last three terms on RHS, the Reynolds stresses and their derivatives, is referred to as turbulence forces. 74 All the terms are divided by the fluid density and are represented in the acceleration unit of 𝑚 𝑠2⁄ . A magnitude comparison has been made between these forces, to gain a deeper understanding of the interplay between the various terms of the radial momentum equation, in different regions of the domain and under different flow conditions. Figure 3.9 shows the magnitude of each force, divided by the fluid density, along the radial direction in the fully developed region of the flow for a stationary pipe (N=0) at different Reynolds numbers. It is evident that the balance between the pressure and turbulence forces dictates the flow field behaviour for the conventional, stationary pipe flow. At all Reynolds numbers, the viscous and convective forces are shown to be inconsequential while the pressure and turbulence forces are the two overwhelmingly dominant forces in the radial momentum equation. Figure 3.9: Accelerations, different forces in equation (2.73) normalized by density, for the stationary wall case (N=0) at different Reynolds numbers. Figure 3.10 compares the different forces at Re= 7000 and x/D = 19 for lower rotation numbers. All the forces were normalized by their respective pipe wall centrifugal force, i.e. rω2. Note that for the sake of visual clarity, the Y-axis limits for N = 0.029 (top left) were chosen to differ from the rest of the graphs. It is observed that at the lowest rotation number (N = 0.029), the pressure and turbulence forces are, by far, the most dominant terms in the radial momentum equation 75 throughout the domain. Therefore, one could reduce the radial momentum equation, i.e. equation (2.74), for such low rotational speeds at Re = 7000, to a balance between the pressure gradient and Reynolds stresses terms with little or negligible error. Doubling the rotational speed to N = 0.057 increases the convective forces especially in regions close to the wall. This is expected as increasing the wall rotational speed extends the reach of the centrifugal force field. However, this increase does not affect the strong dominance of the pressure and turbulence forces within the majority of the domain. It is only for 0.8 < r/R < 0.9 where the convective forces become comparable (and even equal at r/R=0.9) to their turbulence and pressure counterparts. However, the turbulence and pressure forces regain their dominance, owing to their exponential rise, in the vicinity of the wall at radii greater than 0.9, despite the continuous increase of the convective acceleration as the wall is approached. Further increase in the rotational speed (N = 0.14) extends the domain within which the convective forces are comparable to their turbulence and pressure counterparts. It can be seen that for 0.5 < r/R < 0.9 these three forces are comparable. At N=0.29 the convective forces outweigh the turbulence forces for the vast majority of the domain except for in the vicinity of the wall, i.e. 0.95 < r/R. Therefore, for the given Reynolds number of 7000, there is a ten-fold increase in the rotational speed required to shift the force balance from the pressure-turbulence pair at N = 0.029 to the pressure-convective pair at N = 0.29 (excluding the near-wall region), and consequently stabilizing the turbulent flow field. 76 Figure 3.10: Non-dimensionalized accelerations, different forces in equation (2.73) normalized by the rotating wall centrifugal force, for low to intermediate rotation numbers at Re = 7 000. The force balance analysis at higher rotation numbers is illustrated in Figure 3.11. It is seen that the convective forces continue to rise with increasing pipe rotational speed. At N = 0.57 the convective forces match the pressure force, outgrowing the turbulence forces in the process, for almost the entire domain, expect for a close proximity of the wall. At a rotation number of N = 0.86 the convective forces, for the first time, surpass their turbulence counterpart throughout the domain including at the wall and its proximity. This trend intensified with further increases in the rotational speed. At N = 1.43, it can be observed that the convective forces completely outweigh the turbulence forces and counterbalance the pressure force almost single-handedly. The turbulence forces have been diminished entirely by the stabilizing swirling effects, generated by the pipe wall rotation. Similar profiles have been illustrated in Figure 3.12 for a fixed rotation number and various Reynolds numbers. It can be seen that the normalized forces exhibit the same trend irrespective of Reynolds number. This is consistent with the observations made for the stationary pipe in Figure 3.9. 77 Figure 3.11: Non-dimensionalized accelerations, different forces in equation (2.73) normalized by the rotating wall centrifugal force, for intermediate to high rotation at Re = 7000. Figure 3.12: All forces profiles, normalized by their respective pipe wall centrifugal force, for different Reynolds numbers at N=0.29. 78 3.8 An Integral Measure for Evaluating Convective forces vs Turbulence forces In this section we present an integral measure for the pressure force, convective forces, and turbulence forces, and introduce a criterion for identifying dominant contributors to the force balance interplay. To gauge the strength of each force, at a plane of interest, in the fully developed region of the flow, the following integral measure is employed, ‖𝐴‖ = (∫ ( 𝐴(𝑟) − 𝐴 ̅)2 𝑑𝑟𝑟=𝑅𝑟=0)1 2⁄ (3.2) where ?̅? is the average of variable A along the radial direction over the target cross section, defined as, ?̅? =1𝑅∫ 𝐴(𝑟)𝑑𝑟𝑟=𝑅𝑟=0 (3.3) This integral measure is calculated for the pressure force, connective forces, and turbulence forces, and are respectively denoted by ‖𝑃‖, ‖𝑐𝑜𝑛𝑣‖ , and ‖𝑡𝑢𝑟𝑏‖ hereafter. Figure 3.13 presents the ratio of ‖𝑐𝑜𝑛𝑣‖‖𝑃‖ and ‖𝑡𝑢𝑟𝑏‖‖𝑃‖ as a function of the wall rotational speed for different Reynolds numbers. The integrated forces ratios exhibit a similar overall trend, although opposite in magnitude. As rotational speed increases, turbulent forces decrease while convective forces increase. These force ratios become equal at some point. This point, at which the normalized integral measures for the convective and turbulence forces meet, represents a critical rotational speed referred to as the ‘transition point’ in this study. At rotational speeds lower than this critical transition point, the integral measure for turbulence forces is greater than that of the convective forces. On the contrary, at rotational speeds higher than the transition point, the convective forces integral measure surpasses their turbulence forces counterpart. Moreover, the transition rotational speed increases with increasing Reynolds number. This is expected because at higher Reynolds numbers, and correspondingly higher axial flow strength, a stronger swirling flow is needed to suppress a stronger turbulence-driven instability and to stabilize the flow. For instance, the transition rotational speed increases from ω = 0.3 rad/s at Re = 14000 to ω =0.6 rad/s at Re = 28000. A 79 similar trend exists at other Reynolds numbers; there is a linear proportionality between the flow Reynolds number and its transition rotational speed. Interestingly, the graphs of Figure 3.13 collapsed onto a single, universal curve when the wall rotational speed was non-dimensionalized by the axial velocity time scale (𝑅𝜔 𝑈𝑚⁄ ). In this dimensionless space, as shown in Figure 3.14, the transition point between the forces occurs at a unique common rotation number (transition rotation number) of Ntr ≈ 0.45 for all cases tested. This is consistent with the linear proportionality observed between the dimensional transition rotational speed and the flow Reynolds number. Figure 3.13: Non-dimensionalized integral measures of the convective and turbulence forces, normalized by the integral measure for the pressure force, as a function of wall rotational speed for different Reynolds numbers. 80 Figure 3.14: Non-dimensionalized integral measures of the convective and turbulence forces, normalized by the integral measure of the pressure force, as a function of rotation number (N=Rω/Um) for different Reynolds numbers. The transition point between the forces ratios occurs at the same rotation number of Ntr ≈ 0.45. 3.9 Guidelines on Choice of Turbulence Model for Swirling Flows It should be emphasized that the above integral measure analysis was performed for the radial momentum force balance only, indicating a decline in the significance of turbulence in the radial direction with increasing rotation number. This is consistent with the observations made in Figure 3.8; the Reynolds shear stress, −𝜌𝑢𝑧′𝑢𝑟′̅̅ ̅̅ ̅̅ , dropped considerably as the rotation number increased. The focus was put on the radial direction due to its critical role in the particle separation mechanism in swirling separation equipment, such as the hydrocyclone (Averous and Fuentes, 1997; Dai et al., 1999). A comparative study between the k-ε model and RSM is presented below to highlight the importance of the choice of turbulence model, and to propose a criterion for choosing a turbulence model, to ensure an optimum choice in terms of both accuracy and computational cost. Eight further simulation cases were run using the k-ε model at the same rotational speeds as the RSM cases. The k-ε model axial velocity profiles at x/D ≈ 20 are compared against those of the RSM in Figure 3.15. 81 The discrepancy between the axial velocity profiles is insignificant at lower rotation numbers; the maximum difference being less than 6% and 8% for N=0 and 0.14, respectively. Interestingly, at N = 0.57 (> Ntr) the k-ε model predictions start to significantly deviate from those of RSM both in magnitude and overall trend. The maximum difference rises to 35%, while the discrepancy grows in the central core and extends to the near-wall region. Further increase in rotation number leads to a more significant deviation of the k-ε results from those of the RSM. At N = 2.86 the overall trend of the k-ε model axial velocity is completely different from that of RSM (and those reported in literature), see Figure 3.4. Hence, one may conclude that the optimum choice of turbulence model for swirling flows with curved streamlines, becomes imperative when a measure of dimensionless rotational speed exceeds a critical threshold, i.e. the flow rotation number is greater than the transition rotation number, N > Ntr. Figure 3.15: Comparison between the axial velocity predictions of k-ε model and RSM at Re = 7000 for a range of rotation numbers, 0 ≤ N ≤ 2.86. In other words, at rotation numbers lower than Ntr, simpler turbulence models, such as k-ε, are capable of providing satisfactory engineering accuracy; thus avoiding high computational costs of more complex turbulence models, without compromising accuracy. 82 3.10 Summary The effect of swirl on turbulence has been intensively investigated, experimentally and numerically, over the years. The laminarization of turbulent velocity profiles in a rotating pipe has been well established in the literature. Elaborate turbulence models have been shown to provide most accurate descriptions of the classic rotating pipe problem. However, they are very computationally intensive and not affordable to design engineers. Industry will significantly benefit from a cheaper model to replace expensive LES or DNS techniques without compromising accuracy. This work proved that under certain conditions less intensive and more affordable turbulence models could replace their complex and elaborate counterparts, yet providing reasonable accuracy. A criterion was introduced to help determine whether or not a cheaper turbulence model would be sufficient and could be used to predict the flow field accurately. CFD simulations were carried out to give further insight into the effect of swirl on turbulence. The laminarization phenomenon in a turbulent rotating pipe was verified for a range of Reynolds numbers (mean axial velocities) and a variety of rotation numbers (pipe wall rotational speeds). The radial momentum equation was closely analyzed and the evolution of interplay between forces across the pipe was presented as a function of rotation number, for different Reynolds numbers. A new measure, not previously reported, was introduced for analyzing the flow. The presented novel integral measure evaluates the interplay between the pressure, convective, and turbulence forces, and quantifies their contributions to the momentum exchange. It was revealed that there exists a transition rotational speed beyond which the dominance of turbulence over convective forces is reversed. The transition rotational speed was found to be linearly proportional to Reynolds number. A universal transition rotation number of Ntr ≈ 0.45, irrespective of Reynolds number, was found to be the critical point, greater than which the flow transition occurs from pressure-turbulence-governed to pressure-convective-governed. The simulation results provide insight into the significance of different forces, namely the convective, pressure, turbulence and viscous forces, based on the regime of the flow. At rotation numbers greater than the transition Ntr ≈ 0.45, the convective forces outgrew the turbulence and emerged as the dominant competitor for the pressure force. 83 Moreover, the present study makes a further contribution by using the integral measure tool for analyzing the influence of rotation on turbulence. The use of the integral measure led to defining Ntr, which could potentially serve as an indicator of the extent of influence of the choice turbulence model on the accuracy of the CFD predictions. The choice of turbulence model was shown to be crucial for the accuracy of CFD calculations particularly when Ntr is exceeded. While two-equation turbulence models, such as k-ε model, were capable of providing sufficient accuracy for flows with N < Ntr, more complex turbulence models, e.g. RSM, were required when N > Ntr. Guidelines for optimum choice of turbulence model, on the basis of rotation number, were provided. The following observations were made, At very low rotational speeds (N < 0.1), irrespective of Reynolds number, the balance of the radial momentum was governed by the pressure and turbulence forces with little influence from the convective forces. At intermediate rotational speeds (0.1< N < 0.85), irrespective of Reynolds number, the growth of the convective forces enabled them to compete with the turbulence effects in counterbalancing the pressure force. At high rotational speeds (N > 0.85), irrespective of Reynolds number, the convective forces replaced the turbulence forces as the main competitor for the pressure force. At N >1 the convective forces completely outweighed the turbulence forces throughout the domain and, along with the pressure force, dominated the radial momentum balance. The transition rotation number for the turbulence and convective forces was observed to be linearly proportional to Reynolds number. The integral force ratios, normalized by the integral measure for the pressure force, collapsed onto one universal pair of curves when plotted as a function of rotation number. The transition between the turbulence and convective forces occurred at a dimensionless rotation number of Ntr ≈ 0.45 for all cases, irrespective of Reynolds number. The choice of turbulence model is critical for flows with rotational numbers exceeding the transition rotation number (N > Ntr). Flows with rotation numbers < Ntr can be modeled using computationally less expensive models like k-ε model with reasonable accuracy. 84 Chapter 4: Effects of Geometrical and Flow Parameters It has been well established that the separation mechanism is governed by the complex swirling fluid flow inside the hydrocyclone. As mentioned in Chapter 1, in most studies concerned with geometry modification the emphasis has been put on the particle separation performance (the outcome) rather than the response of the flow field (cause) to the changes. A better understanding of the flow field can greatly facilitate the development of more reliable CFD tools designed to optimize hydrocyclones for specific applications. In this chapter, the hydrodynamic response to a cone angle change is closely analyzed. The pertinent details will be discussed and the physics behind changes in flow field features are explained. Three hydrocyclone geometries with different cone angles and conical body sizes were considered. The hydrocyclones’ dimensions are presented in Table 4.1 where subscripts c, i, v, o, u and con, respectively refer to cylindrical section, inlet, vortex finder, overflow, underflow, and conical section. The dimensions of case A are identical to the hydrocyclone experimentally studied by Mackenzie (2014). The inlet is a circular pipe tangentially connected to the cylindrical body of the hydrocyclone. Figure 4.1 portrays a comparison between the three hydrocyclone geometries, with all the dimensions identical, except for the cone angles and cone lengths. A tetrahedral grid was employed for the computational domain of all cases. Figure 4.2 depicts the grid used for case A. Similar tetrahedral grids were employed for cases B and C. Table 4.1: Hydrocyclones dimensions in mm. Case D di do du Lv Lcyl Lcon θ(°) A 102 19 12.7 6 63.5 102 405. 6.76 B 102 19 12.7 6 63.5 102 272 10 C 102 19 12.7 6 63.5 102 179 15 85 (a) (b) Figure 4.1: (a) The hydrocyclone general geometry; (b) overlay of the three cone angle geometries, cases A (red), B (green), and C (blue). Figure 4.2: Tetrahedral grid of case A. Similar grid was used for cases B and C. 86 4.1 Numerical Scheme and Boundary Conditions Numerical simulations of three-dimensional, turbulent, and incompressible flow of water in the hydrocyclone were performed by means of CFD package ANSYS FLUENT. Both Reynolds Stress and Large Eddy Simulation models were employed for turbulence modelling and their results were compared with the experimental data. The quadratic pressure-strain was chosen for the RSM due to its superior performance compared to LPS, as discussed in Chapter 1. The high-order QUICK scheme was employed for the discretization of the momentum equations. The inlet boundary condition was set to velocity inlet, whereas pressure outlet boundary conditions were assigned to the underflow and overflow exits. The averaged pressure at the overflow was monitored, and adjusted when needed, to ensure desired reject rate (RR) and split ratio in accordance to the experimental data. In order to preserve the swirling nature of the flow, the radial pressure gradient was allowed at the outlet orifices. No air-core formation was observed in the experimental study of the hydrocyclone at RR=50%; thus, no air core was defined in the simulations. The grid used in the RSM simulations contains approximately 350,000 cells, whereas the LES grid consists of nearly 1,200,000 cells. A time step size of ∆𝑡 = 10−5𝑠 was chosen for the LES simulations where the RSM solution was interpolated to be used as the initial solution. For the RSM simulations, the time step size was ∆𝑡 = 10−4𝑠, and k-ε solution was used for initialization. The simulations were run for several residence times to ensure that the flow field was established and no significant changes of monitored flow variables (and their time-averaged values) were observed. 4.2 Experimental Validation The numerical solutions obtained by both RSM and LES were compared against LDV measurements. Figure 4.3 presents the tangential velocity profiles at four different axial heights as a function of radial position normalized by the local hydrocyclone radius (R). The velocity inlet is set to 2.6 m/s (corresponding to a mass flow rate of 0.72 kg/s) and the pressure values at the outlets result in a reject rate (RR) of approximately 50%. It can be observed that the LES and RSM results give similar predictions within the free-vortex region (outer core) and are reasonably close to the LDV measurements. However, the RSM under-predicts the tangential velocity within the forced-87 vortex region (inner core), whereas the LES solutions are in excellent agreement with the experimental data. The under-prediction of the tangential velocity by RSM is consistent by the observations of other authors, discussed in Chapter 1. It should be noted that there is an uncertainty of 4 mm of the radial position for the experimental data. For more detail about the flow loop and error analysis refer to Mackenzie (2014). Figure 4.3: Clockwise from top left; tangential velocity profiles at z1=71.9, z2=91, z3=110.1, and z4=129.1 mm below the inlet vs. the radial position normalized by the local wall radius. 88 Figure 4.4: Clockwise from top left; tangential velocity profiles at z5=148.2, z6=167.3, z7=186.4, and z8=205.5 mm below the inlet vs. the radial position normalized by the local wall radius. A similar comparison for the axial velocity profiles is presented in Figure 4.5. It can be seen that the predictions of both turbulence models are in reasonably good agreement with the experimental data. The LES model provides a closer match both in terms of magnitude and trend particularly within the forced-vortex region. Similar agreement exists for both tangential and axial velocity components at other heights across the hydrocyclone. 89 Figure 4.5: Clockwise from top left; tangential velocity profiles at z=71.9, z=91, z=110.1, and z=129.1 mm below the inlet vs. the radial position normalized by the local wall radius. To further test the simulation results, the area-weighted average of the axial velocity was computed at various planes of different axial heights, and was compared with the analytical solution. From the continuity we know that the averaged axial velocity at an axial plane (perpendicular to the cyclone axis) must satisfy the following equation, 𝑉𝑎𝑥𝑖𝑎𝑙 =?̇?𝜌𝐴 (4.1) where ?̇? is the mass flow rate and A is the cross-sectional area. Note that the cross-sectional area within the conical body is a function of distance from the top of the conical body, Z, i.e. the border of the cylindrical and conical bodies, and the cone angle, θ. The cross-sectional area of a plane can be expressed as 𝐴 = 𝜋𝑟2 with 𝑟 = 𝑅𝑐𝑦𝑙 − 𝑍 𝑡𝑎𝑛𝜃. Figure 4.6 shows that the numerical and analytical solutions are in excellent agreement. This serves as additional validation for the conservation of mass in the numerical calculations. Note that the averaged axial velocity increases by almost one order of magnitude over the length of the conical section. 90 Figure 4.6: Axial velocity; theoretical solutions vs. numerical solutions. Vorticity contours on both an axial plane and a plane perpendicular to the hydrocyclone axis are illustrated in Figure 4.7. It is observed that the vorticity magnitude in the outer core is several orders of magnitude smaller than its values in the inner core. This indicates that the flow is irrotational in the outer core (free vortex), and rotational in the inner core (forced vortex). This will be discussed in more detail in subsections 4.3.4 and 4.4. 91 (a) (b) Figure 4.7: Vorticity magnitude contours at (a) mid-plane x=0, and (b) plane z1 for case A (1/s). 4.3 Effect of Cone Angle 4.3.1 On Tangential Velocity Tangential velocity profiles at different heights are illustrated in Figure 4.8 for all three cases. To ensure geometric similitude, the velocity profiles are extracted at planes with the same cross-sectional areas, rather than the same axial heights, in the conical body. The cross-sectional area of planes at the same axial heights is different owing to the change of cone angle. Thus, the axial heights were chosen in such a way that the ratio of the distance from the top of the conical body to the length of the conical body remains the same across different cases, for instance (h1/Lcon)A = (h1/Lcon)B =(h1/Lcon)c. For example, if h1 is the axial distance between the conical-cylindrical borderline and the line at which the tangential velocity profile is extracted in a case with a cone angle θ1, the corresponding height, h2, in a cone angle θ2 case will have to satisfy the following relationship; ℎ1 tan 𝜃1 = ℎ2 tan 𝜃2. 92 Figure 4.8: Tangential velocity profiles at (a) z=71.9, (b) z=91, (c) z=110.1, and (d) z=129.1 mm below the inlet vs. the radial position normalized by the local wall radius. It can be seen that changes in tangential velocity, and consequently in the centrifugal field, as a result of changes in the cone angle are smallest within the free-vortex region, including the wall proximity. The flow local acceleration in this region is primarily dictated by the near-zero vorticity field shown in Figure 4.7. Consequently, the pressure field passively adjusts itself with the translational acceleration, which could be obtained theoretically from the conservation of mass. Given that the planes at which the profiles were extracted are of the same cross-sectional area, changes in the axial velocity field should be minimal, and only due to the general reduction in the cross sectional area. Hence, the tangential velocity field in the irrotational region of the flow remains relatively insensitive to changes in cone angle. The discrepancy between the tangential velocity profiles becomes more pronounced towards the inner core of the free-vortex region, i.e. towards the hydrocyclone axis. The most significant changes in the tangential velocity field with cone angle occur within the forced-vortex core. The tangential velocity peak increases the most dramatically with an increase in cone angle. The relative difference between the tangential velocity profiles for cases A (θ =6.76°) and C (θ =15°) has been depicted in Figure 4.9. The tangential velocity relative difference (TVRD) is defined as (Vθ, C-V θ, A )/ V θ, A. It is seen that at both heights the TVRD remains well under 20% 93 for the majority of the domain when y/R > 0.7. However, it gradually increases by a small but steady slope as y/R decreases (0.2 < y/R < 0.6). The most dramatic response to the cone angle change is observed in the vicinity of the hydrocyclone axis, including the entire forced-vortex core. The smallest maximum increase occurs inside the cylindrical body at the axial location z1. Figure 4.9: Relative difference in the tangential velocity profiles of the cases C (θ = 15) and A (θ = 6.76) at z1 and z3. 4.3.2 On Pressure Figure 4.10 presents the distribution of the radial pressure drop, p-pinlet, at z3 (110.1 mm below the inlet). The pressure drop within the free-vortex core, especially in the proximity of the wall, is negligible compared to that in the forced-vortex core. This is due to the near-zero vorticity field and small changes in the velocity magnitude within the free-vortex core, as presented in section 4.3.1. The pressure drop grows sharply and becomes significant within the forced-vortex core, reaching its maximum at the hydrocyclone axis in all three cone angle cases. The same trend was observed at other heights. The interplay between the pressure force and the acceleration field is more involved in this region of the flow as the vorticity field plays a critical role in the flow dynamics. 94 It should be noted that the residence time of particles at the free-vortex and forced-vortex regions of the flow are considerably different. On their downward motion towards the underflow, the fluid particles initially experience the vorticity-free region, where pressure changes with respect to the inlet pressure are negligible. Later, a significant portion of the particles would experience the high vorticity regions in their motion towards either the overflow or underflow. The sequence of events could be described as follows. The flow is initially driven by the pressure force inside the hydrocyclone, where the tangential inlet converts the translational motion to the rotational motion. As the rotational flow moves further down in the hydrocyclone, the fluid particles angular velocity, and consequently their vorticity, increases due to the reduction in the cross-sectional area in the conical section. The multiplication of the tangential velocity and the radius of rotation should remain nearly constant (see section 4.4.1 for further analysis). Eventually, the fluid elements experience the high rotational flow prior to reporting to either the underflow or the overflow. The high swirl in the flow, followed by an increase in tangential velocity, results in a higher pressure drop in the forced-vortex region of the flow, as shown in Figure 4.10. Figure 4.10: Profiles of the pressure drop at z3. Radial position is normalized by the local wall radius. While the pressure drop remains relatively constant outside the forced-vortex core (r/R > 0.2), the pressure drop profiles start to deviate from one another inside the forced-vortex core, i.e. at 95 dimensionless radial positions of r/R < 0.2. The discrepancy peaks at the centre where the pressure drop grows in magnitude for case C (Δp > 200 kPa), more than doubles its counterpart in case A. This observation is consistent with the tangential velocity profiles presented in section 4.3.1. Figure 4.11: Pressure drop at z3 normalized by the pressure difference between underflow and overflow. Figure 4.11 presents a dimensionless form of the pressure drop of Figure 4.10, non-dimensionalized by the pressure difference between the underflow (pu) and overflow (po). Figure 4.12 depicts a new dimensionless account of the pressure field. The dimensionless pressure drop profiles of Figure 4.11 have been normalized by the squared ratio of the inlet velocity to the velocity magnitude at the plane of interest. The commonly used scaling relationship is given as 𝐸𝑢 = 𝑓𝑢𝑛𝑐(𝑅𝑒), where Eu is the Euler numer defined as 𝐸𝑢 = ∆𝑝 𝜌𝑉2⁄ , and 𝑅𝑒 = 𝜌𝑉𝑟𝑒𝑓𝐷 𝜇⁄ , where D is the hydrocyclone characteristic size and 𝑉𝑟𝑒𝑓 is the characteristic velocity of the flow (inlet velocity). The new dimensionles pressure drop can be interpreted as the ratio of two different Euler numbers; a local Euler number defined as (p-pinlet) / V2local and a global Euler number defined as (pu-po) / V2inlet. This new normalization causes the dimensionless pressure drop profiles to almost perfectly collapse onto a single graph, irrespective of the cone angle magnitude. This suggests a universal normalized dimensionless pressure drop could be introduced for the 96 hydrocyclone, regardless of cone angle value. This new dimensionless group could potentially be useful in scaling dissimilar hydrocyclones. Figure 4.12: The distribution of the ratio of the local Euler number to the global Euler number at z3. Figure 4.13: Pressure difference ratio (PDR) for the three cone angle cases. 97 The ratio of the overflow pressure drop to the underflow pressure drop (relative to the inlet pressure) is called the Pressure Difference Ratio (PDR). The PDR is an important operational parameter for the hydrocyclone classification and is defined as, 𝑃𝐷𝑅 =𝑝𝑖𝑛𝑙𝑒𝑡 − 𝑝𝑜𝑝𝑖𝑛𝑙𝑒𝑡 − 𝑝𝑢 (4.2) Figure 4.13 shows the changes of PDR with cone angle. It is seen that an increase in cone angle leads to a higher PDR. However, changes in PDR with cone angle appear to diminishes with further increases in cone angle; a similar trend was reported by Saidi et al., (2013). An increase in cone angle is followed by a decrease in the conical length, owing to the unchanged underflow opening size. This results in a shorter distance for fluid particles to travel to the underflow. Therefore, the rate of change of pressure drop between the inlet and the underflow increases mildly with an increase in cone angle. On the other hand, given that the reject rate is kept constant among the cases, a higher pressure drop is required between the inlet and the overflow to satisfy the desired flow split; i.e. the rate of change of pressure drop increases sharply with an increase in cone angle. Hence, the increase in the cone angle tends to increase the numerator faster than then the denominator in equation (4.2), leading to an increase in the PDR which the data suggest eventually levels off. 4.3.3 On Viscous Dissipation The viscous dissipation function was computed for each case in order to analyze the effect of cone angle on energy loss. The viscous dissipation, Φ, for an incompressible flow of a Newtonian fluid is defined as, Φ = 2𝜇 [(𝜕𝑢𝜕𝑥)2+ (𝜕𝑣𝜕𝑦)2+ (𝜕𝑤𝜕𝑧)2]+ 𝜇 [(𝜕𝑣𝜕𝑥+𝜕𝑢𝜕𝑦)2+ (𝜕𝑤𝜕𝑦+𝜕𝑣𝜕𝑧)2+ (𝜕𝑢𝜕𝑧+𝜕𝑤𝜕𝑥)2] (4.3) where u, v, and w, are velocity components in the x, y, and z directions, respectively, and μ is the molecular viscosity. Volume-averaged values of the viscous dissipation function were computed 98 for the entire hydrocyclone, ΦC, as well as for a central core of the cyclone, ΦUT, formed by a cylinder with a radius of r=3mm (same size as the underflow opening) extending from the underflow to the overflow. The dissipation function values for the two volumes, as well as the ratio of the volume of the underflow tube to the entire volume of its respective cyclone, are given in Table 4.2. Table 4.2: Volume-averaged values of the dissipation function (ΦUT for r=3 mm) Case ΦUT (underflow tube) ΦC (entire cyclone) ΦUT / ΦC VUT/VC θ=6.6 (reference) 5.37 × 102 1.37 × 104 25.63 8.03 × 10−3 θ=10 7.458 × 102 2.4552 × 104 32.92 7.86 × 10−3 θ =15 1.14 × 103 3.92 × 104 34.47 7.48 × 10−3 It is observed that while in all cases the underflow tube forms a very small fraction of the entire hydrocyclone volume (only about 2%), the volume-averaged viscous dissipation in the tube is significantly greater than its counterpart in the entire cyclone. This indicates that the vast majority of energy loss occurs within the central core (forced-vortex core), consistent with observations made for the tangential velocity and pressure profiles. Moreover, the largest viscous dissipation occurs in the θ=15 case which has the smallest volume among the three cases. This may sound contradictory as one would expect that a decrease in volume, accompanied by a decrease in the wall surface area, should result in reduced friction and consequently less viscous dissipation. However, it should be noted that in addition to the surface area, the kinetic energy is another factor which plays an influential role in energy loss due to the viscous dissipation. The tangential velocity component, as the dominant term in the kinetic energy, increases significantly with cone angle, as observed in Figure 4.8. This increased kinetic energy overrides the effect of the surface area reduction and ultimately results in a higher dissipation than that in cases with a smaller cone angle, i.e. larger surface area and smaller kinetic energy. 4.3.4 On Vorticity The vorticity magnitude profiles extracted along the radial direction at an axial location z3 are presented in Figure 4.14. For all three cases, the vorticity magnitude is negligible for the majority 99 of the radial domain (r/R > 0.2), the free-vortex region, except for the regions adjacent to the wall, where the flow becomes rotational due to the viscous boundary layer. Note that the negative values for the near-wall region, where the flow becomes rotational, is due to increased shear strain rate in close proximity to the wall. Similar observations were made by other authors (Dyakowski and Williams, 1993; Fisher and Flack, 2002). Moreover, the vorticity magnitude starts to grow at approximately r/R ≈ 0.2 and increases dramatically within the forced-vortex core (r/R < 0.2). The vorticity reaches its peak value in the close proximity of the hydrocyclone axis. Additionally, the higher the cone angle, the higher the vorticity magnitude in the inner core of the hydrocyclone. Figure 4.14: Vorticity magnitude profiles at z3 vs. dimensionless radial position. It can be seen that the vorticity field is near-zero in the majority of the outer region of the flow field, i.e. r/R > 0.2, corresponding to the free-vortex region of the tangential velocity field. The imperfect irrotationality is due to the presence of viscous forces in that region which lead to a slight deviation of the free-vortex region from the ideal potential flow. This is in line with the observations reported by Fisher and Flack (2002). Note that the negative values for the vorticity in the near-wall region where the flow becomes rotational is due to an increased shear strain rate in the close proximity of the wall. Similar observations have been made by other authors (Fisher and Flack, 2002; Dyakowski and Williams, 1993). At a radial position around r = 10 mm, the vorticity magnitude starts to grow and undergoes a transition from the free-vortex type of flow to 100 a forced-vortex type of flow. Further towards the hydrocyclone axis, the vorticity magnitude reaches its maximum value inside the forced-vortex region, where the flow becomes fully rotational. It can be seen that there is a direct relationship between the cone angle and the vorticity magnitude; the higher the cone angle, the stronger the vorticity field, and consequently the greater the tangential velocity. A similar trend exists for other axial locations, not shown here. 4.3.5 On G-Force The tangential velocity field determines the distribution of the centrifugal field acting on the disperse phase. The equation of motion for a spherical solid particle in the radial direction (assuming the slip between the fluid elements and solid particle in the 𝜃-direction is negligible) takes the following form, (𝜌𝑝 − 𝜌𝐹)∀𝑉𝜃2𝑟− 𝐶𝑑𝐴 (12𝜌𝐹𝑉𝑟,𝑠𝑙𝑖𝑝2 ) = ∀𝜌𝑝𝑑𝑉𝑟,𝑠𝑙𝑖𝑝𝑑𝑡 (4.4) where ∀=𝜋6𝑑𝑝3 is the volume of a spherical particle with diameter 𝑑𝑝, 𝐴 =𝜋4𝑑𝑝2 is the cross-sectional of the particle, 𝑉𝑟,𝑠𝑙𝑖𝑝 is the slip velocity between the fluid elements and the solid particle in the radial direction, 𝐶𝑑 is the drag coefficient based on the slip velocity, and 𝜌𝑝 , 𝜌𝐹 are respectively, the density of the fluid and density of the particle. Rearranging of equation (4.4) leads to, Λ𝑉𝜃2𝑟− 𝑉𝑟,𝑠𝑙𝑖𝑝2 = Λ𝜌𝑝(𝜌𝑝 − 𝜌𝐹)𝑑𝑉𝑟,𝑠𝑙𝑖𝑝𝑑𝑡 (4.5) where Λ =43(𝜌𝑝−𝜌𝐹)𝜌𝐹𝑑𝑝𝐶𝑑, which can be viewed as a characteristic length scale on the particle size. Given a small relative velocity between the fluid elements and solid particle, the drag coefficient is 𝐶𝑑 = 𝐶𝑑(𝑅𝑒𝑟𝑒𝑙) = 24 𝑅𝑒𝑟𝑒𝑙⁄ , where 𝑅𝑒𝑟𝑒𝑙 = 𝜌𝐹𝑉𝑟,𝑟𝑒𝑓𝑑𝑝 𝜇⁄ (𝑉𝑟,𝑟𝑒𝑓 is the characteristic relative velocity between particle and fluid elements), leading to, 101 Λ =43(𝜌𝑝 − 𝜌𝐹)𝜌𝐹𝑑𝑝𝐶𝑑= [(𝜌𝑝 − 𝜌𝐹)𝑑𝑝218𝜇]⏟ 𝜏𝑝𝑉𝑟,𝑟𝑒𝑓 (4.6) The terms in the brackets represents the relaxation time of the solid particles in the flow. Rendering equation (4.5) dimensionless by the flow characteristic length scale D, characteristic velocity 𝑉𝑟,𝑟𝑒𝑓 and characteristic time scale 𝜏𝐹 = 𝐷 𝑉𝑟,𝑟𝑒𝑓⁄ , we obtain, 𝜏𝑝𝜏𝐹?̃?𝜃2?̃?− ?̃?𝑟,𝑠𝑙𝑖𝑝2 =𝜏𝑝𝜏𝐹𝜌𝑝(𝜌𝑝 − 𝜌𝐹)𝑑?̃?𝑟,𝑠𝑙𝑖𝑝𝑑?̃? (4.7) where the tilde notation denotes dimensionless quantities. The ratio of the particle time scale and the flow time scale is commonly known as the Stokes number. This equation shows that the centrifugal field strength is weighted by the Stokes number. The radial distributions of the centrifugal field for the three cone angle cases have been illustrated in Figure 4.15. It is clear that the magnitude of the centrifugal field is affected the least by the cone angle changes in the free-vortex region. As the hydrocyclone axis is approached, at radii r/R < 0.3 the graphs start to deviate. The deviation grows sharply with increases in the radial position and peaks inside the forced-vortex core at r/R ≈ 0.1. The centrifugal field undergoes a maximum change of a five-fold increase from case A to case C. This behavior is in line with the observations made for the tangential velocity and pressure field in sections 4.3.1, 4.3.2, and 4.3.4. 102 Figure 4.15: Radial distribution of centrifugal field for the three different cone angle cases at z3. 4.4 Force Balance Analysis To gain further insight into the underlying force balance in the flow, we closely examine the conservation of momentum equation, the Navier-Stokes equation, in the following form (assuming that gravity acts in the –z direction), ∇⃗ (𝑝𝜌+?⃗? . ?⃗? 2+ 𝑔𝑧⏟ 𝐵) = ?⃗? × Ω⃗⃗ −𝜇𝜌∇⃗ × ?⃗? (4.8) where p is the instantaneous pressure, ?⃗? is the instantaneous velocity field and ?⃗? = ∇⃗ × ?⃗? is the instantaneous vorticity vector. The term within the parentheses, B, is the Bernoulli function, which represents the total energy per unit mass of the fluid. The first term on the RHS represents rotationality effects and is referred to as the vortex force, whereas the second term on the RHS represents the viscous forces. It is worth remembering that under irrotationality conditions, i.e. in the free-vortex region, B remains constant throughout the irrotational flow domain. Employing the Reynolds decomposition, the averaged momentum equation takes the form, 103 ∇⃗ (𝑃𝜌+?⃗? . ?⃗? 2+ 𝑔𝑧⏟ ?̅?) = ?⃗? × Ω⃗⃗ −𝜇𝜌∇⃗ × Ω⃗⃗ − ∇⃗ . 𝜏𝑖𝑗𝑡̿̿ ̿ (4.9) where ?⃗? and Ω⃗⃗ refer to the mean velocity and vorticity fields respectively, 𝑃 denotes the mean pressure field, and 𝜏𝑖𝑗𝑡̿̿ ̿ is the Reynolds stress tensor given by 𝜏𝑖𝑗𝑡̿̿ ̿ = −𝜌𝑢𝑖′𝑢𝑗′̅̅ ̅̅ ̅̅ . The term within the parentheses in equation (4.9) is called the averaged Bernoulli function, ?̅?. The radial components of equation (4.9) along the radial direction at an axial location z3 are shown in Figure 4.16. It is observed that the viscous forces and the turbulence contribution to the total force balance are negligible. However, the turbulence force is observed to be several orders of magnitude greater than the viscous forces throughout the domain, with the forced-vortex core experiencing the highest turbulence effects. Therefore, the hydrocyclone flow field could be interpreted as being composed of two regions: (i) a free-vortex region where ∇⃗ ?̅? = 0, and (ii) a forced-vortex region where ∇⃗ ?̅? = V⃗ × Ω⃗⃗ . The velocity and vorticity vectors are expanded, for further analysis, in the following section. Figure 4.16: (right) Bernoulli gradient, vortex force, viscous force, and turbulence; (left) viscous and turbulence forces, at z3 for case A. The forces represent the terms in equation (4.9) multiplied by density, ρ=998.2 kg/m3. 104 4.4.1 Vortex Force The r-component of the vortex force in the cylindrical coordinate can be written as, (?⃗? × Ω⃗⃗ )𝑟= 𝑉𝜃Ω𝑧 − 𝑉𝑧Ω𝜃 (4.10) where ?⃗? (𝑟, 𝜃, 𝑧) = 𝑉𝑟(𝑟, 𝜃, 𝑧)?̂?𝑟 + 𝑉𝜃(𝑟, 𝜃, 𝑧)?̂?𝜃 + 𝑉𝑧(𝑟, 𝜃, 𝑧)?̂?𝑧 is the mean velocity field and Ω⃗⃗ (𝑟, 𝜃, 𝑧) = Ω𝑟(𝑟, 𝜃, 𝑧)?̂?𝑟 + Ω𝜃(𝑟, 𝜃, 𝑧)?̂?𝜃 + Ω𝑧(𝑟, 𝜃, 𝑧)?̂?𝑧 is the mean vorticity field. Figure 4.17 depicts a comparison between the magnitudes of the two components of the radial vortex force. The first term, 𝑉𝜃Ω𝑧, is shown to be the significant contributor to the vortex force. The overwhelming dominance of VθΩz is physical as Vθ (>Vz) and Ωz (>Ωθ) are by far the most dominant components of the velocity and vorticity fields. Consequently, their multiplication is several orders of magnitudes greater than UzΩθ. Therefore, one could reduce the radial vortex force to (?⃗? × Ω⃗⃗ )𝑟= 𝑉𝜃Ω𝑧. Figure 4.17: A comparison between the two components of the vortex force at z3 for case A. A closer examination of the vortex force can be performed by further breaking down its only influential term, Vθ Ωz, The axial component of the vorticity field is expressed in the cylindrical coordinates as below, 105 Ω𝑧 =𝜕𝑉𝜃𝜕𝑟+𝑉𝜃𝑟−1𝑟𝜕𝑉𝑟𝜕𝜃 (4.11) Due to the weak circumferential asymmetry, it is expected that the last term on the RHS is much smaller than the other two terms, as shown in Figure 4.18. Figure 4.18: Comparison between the axial vorticity terms of equation (4.11) at z3 for case A. The three terms of the axial vorticity, i.e. the terms on the RHS of equation (4.11), were calculated separately, and are shown in Figure 4.18. The first two terms, ∂Vθ/∂r and ∂Vr/∂θ, are comparable and significantly greater than (1/r) ∂Vr/∂θ. Hence, Ω𝑧 ≈𝜕𝑉𝜃𝜕𝑟+𝑉𝜃𝑟 (4.12) The radial vortex force is accordingly reduced to the following, (?⃗? × Ω⃗⃗ )𝑟≈ 𝑉𝜃 (𝜕𝑉𝜃𝜕𝑟+𝑉𝜃𝑟) (4.13) 106 4.4.2 Bernoulli Gradient The Bernoulli function in the radial direction is composed of the sum of the pressure (P/ρ) and the kinetic energy (?⃗? .?⃗? 2). Hence, its gradient is expressed as, 𝜕𝐵𝜕𝑟=𝜕𝜕𝑟(𝑃𝜌) +𝜕𝜕𝑟(?⃗? . ?⃗? 2) (4.14) The gradient of the kinetic energy can be expanded to, 𝜕𝜕𝑟(?⃗? . ?⃗? 2) = 𝑉𝑟𝜕𝑉𝑟𝜕𝑟+ 𝑉𝜃𝜕𝑉𝜃𝜕𝑟+ 𝑉𝑧𝜕𝑉𝑧𝜕𝑟 (4.15) The second term on the RHS of equation (4.15) is expected to be the most important term due to the presence of the tangential velocity and its radial derivative. Figure 4.19 shows the magnitude of the three terms in equation (4.15). The overwhelming dominance of the circumferential term, Vθ ∂Vθ/∂r, is clearly observed. Thus, the kinetic energy part of the Bernoulli gradient, equation (4.15), can be reduced to, 𝜕𝜕𝑟(?⃗? . ?⃗? 2) ≈ 𝑉𝜃𝜕𝑉𝜃𝜕𝑟 (4.16) Substituting equations (4.13) - (4.16) into equation (4.9) gives, 𝜕𝜕𝑟(𝑃𝜌) + 𝑉𝜃𝜕𝑉𝜃𝜕𝑟= 𝑉𝜃 (𝜕𝑉𝜃𝜕𝑟+𝑉𝜃𝑟) (1) The terms containing the radial derivative of the tangential velocity, Vθ ∂Vθ/∂r, cancel each other out from the LHS and RHS, reducing the radial momentum balance to the following, 𝜕𝜕𝑟(𝑃𝜌) = 𝑉𝜃𝑉𝜃𝑟 (4.17) which is the well-known balance between the radial pressure gradient and the centrifugal force, ∂P/∂r=ρ Vθ2/r. 107 Finally, the radial gradient of Bernoulli has been plotted in Figure 4.20 for the three cone angle cases at the axial position z3. It is shown the Bernoulli gradient magnitude goes up sharply with cone angle. The gradient peak increases by almost 100% from case A to case B, and similarly doubles from case B to case C. This trend is attributed to the rapid changes in the pressure field and dramatic increase of the radial pressure gradient discussed in section 4.3.2. Figure 4.19: The radial distribution of the kinetic energy gradient terms, LHS of equation (4.15), at z3 for case A. 108 Figure 4.20: The radial distribution of Bernoulli gradient at z3 for all cone angle cases. 4.5 Summary CFD simulations of the flow field of three hydrocyclones were performed using both RSM and LES turbulence models. The effects of changes in the design variables on the flow field was investigated. Particular attention was paid to cone angle due to its imperative role in particle separation. The effect of cone angle on the flow characteristics and hydrodynamics was closely analyzed. The following observations were made: (i) The velocity field was compared with LDV experimental data. It was found that both the LES and RSM predictions for the tangential velocity within the irrotational, free-vortex core were accurate. However, in the forced-cortex core, LES was shown to produce a superior performance, whereas the RSM under-predicted the tangential velocity. (ii) Increasing cone angle resulted in a significant increase in the swirl velocity in the forced-vortex core, where sharpest increases occurred in the maximum swirl velocity. However, the free-vortex region exhibited less sensitivity to changes in cone angle magnitude. The same statement applies to the vorticity magnitude in the domain. 109 (iii) A new dimensionless pressure drop was introduced causing the pressure profiles of different cases to collapse onto a single curve. This normalized pressure group can prove helpful in scaling hydrocyclones. (iv) The viscous dissipation function and energy loss were shown to be the highest for the largest cone angle geometry. This was surprising as the largest cone angle geometry resulted in the lowest total volume, smallest wall surface area and consequently lowest friction. However, the significant increase in the tangential velocity, and thus in the velocity magnitude, in the largest cone angle case proved to be a more decisive factor than the surface area. This led to the largest energy loss among cases overriding the surface area effect in the process. (v) The Navier-Stokes equations were presented in the form of a balance between the Bernoulli function and the vorticity effects. The contribution of each term to the total force balance was investigated numerically. It was revealed that the radial momentum balance was governed by a balance between the Bernoulli gradient and the vortex force, ∇⃗ (?̅?) = ?⃗? ×Ω⃗⃗ . (vi) It was shown that the radial vortex force was overwhelmingly dependent on the product of the axial vorticity and the tangential velocity, (?⃗? × Ω⃗⃗ )𝑟≈ 𝑉𝜃 (𝜕𝑉𝜃𝜕𝑟+𝑉𝜃𝑟). (vii) The radial Bernoulli gradient was shown to be mainly dependent on radial pressure gradient and the product of the tangential velocity and its radial derivative, 𝜕𝐵𝜕𝑟=𝜕𝜕𝑟(𝑃𝜌) + 𝑉𝜃 (𝜕𝑉𝜃𝜕𝑟+𝑉𝜃𝑟). 110 Chapter 5: Effects of Vortex Finder Diameter, Number of Inlets and Inlet Positioning The effects of vortex finder diameter and the number of inlets on the hydrocyclone flow field is investigated and presented in this chapter. As discussed in Chapter 1, the effects of the vortex finder and the number of inlets on the separation performance has been studied by several authors. The focus in this study is on the influence of these design variables on the force balance and velocity field specifically, which can be linked to the separation performance. A force balance analysis is performed to shed light on the kinetic response of the flow field to changes in these design variables. 5.1 Effects of Vortex Finder Diameter 5.1.1 On Tangential Velocity Table 5.1 presents the dimensions of different geometries used for this part of the study. Tangential velocity profiles for the vortex finder diameter cases (cases A, B, C) have been illustrated in Figure 5.1. It is observed that the tangential velocity is insensitive with changes in the vortex finder for r/R > 0.5, where the flow is mainly irrotational. The tangential velocity profiles begin to deviate from one another at radii farther away from the wall, i.e. r/R ≈ 0.5. The discrepancy increases considerably at smaller radii and specifically within the forced-vortex core, i.e. r/R < 0.2. It can be seen that the tangential velocity peak is inversely proportional to the ratio of the size of the overflow (vortex finder) to the underflow openings, do/du. Decreasing do/du from 4 in case C to 1 in case B results in an approximately 100% increase in the tangential velocity peak at the presented axial heights. Note that all the flow and geometrical parameters are kept constant, expect for the overflow diameter. Given the fact that the mass influx through the overflow exit is constant across the cases, a smaller vortex finder diameter translates into a higher axial velocity field to maintain the same flow rate. The effect of this increase is most significant on fluid elements near and under the vortex finder, and diminishes as the fluid particles travel away from the vortex finder. The elevated axial 111 velocity field brings about an increase in the tangential velocity field. The physics behind the changes in the tangential velocity will be explained in more detail in section 5.1.2. Table 5.1: Dimensions of different geometries in mm. Case D di do du Lv Lcyl Lcon θ(°) A 102 19 6.35 6 63.5 102 405. 6.76 B 102 19 12.7 6 63.5 102 405 6.76 C 102 19 25.4 6 63.5 102 405 6.76 Figure 5.1: Tangential velocity profiles for the three vortex finder diameter cases at axial distances of (top left) z1=71.9 mm, (top right) z2=91.0 mm, (bottom left) z3=110.1 mm, and (bottom right) z4=129.1 mm below inlet. Radial position is normalized by the respective local wall radius. 5.1.2 On Pressure Figure 5.2 illustrates the change of the PDR, defined by equation (4.2), with changes in the ratio of the overflow and underflow diameters. It can be observed that the PDR declines with increases 112 in the ratio. A sharp drop is seen from do/du=1 to do/du=2, whereas the rate of change from do/du=2 to do/du=4 is mild. This suggests that a geometry with higher ratio of the overflow to underflow diameter would be more suitable for fractionation of smaller coarser particles. Figure 5.2: Pressure Difference Ratio vs. the ratio of overflow diameter to underflow diameter for cases A, B, and C. 5.1.3 On Vorticity A comparison between the axial vorticity profiles of case A - C at different heights are presented in Figure 5.3 -Figure 5.5). It can be seen that the axial vorticity magnitude in case B is drastically larger, at all heights, than its counterparts in the other two cases. However, the axial vorticity magnitude dominance of case B over that of the benchmark geometry, i.e. case A, shrinks from five-folds at z1 (inside the cylindrical body) to three-folds at z10 (bottom of the conical section). This is partly due to the increase of the vorticity magnitude within the conical section. Narrowing of the cross-sectional area in the conical section leads to greater rotation, and consequently to higher vorticity. More importantly, however, the larger distance from the vortex finder is responsible for the shrinking gap between the peaks from inside the cylindrical body to the top of the conical body. At heights near the vortex finder, i.e. inside the cylindrical body mostly and the top part of the conical body, the effects of the vortex finder diameter prevails. The farther from the 113 vortex finder the smaller its effect on the vorticity and flow field. Therefore, the discrepancy between the peak vorticity magnitudes in the three cases decreases at locations far away from the vortex finder, as the effect of the size of the overflow opening starts to diminish. Moreover, it is seen in Figure 5.3 -Figure 5.5) that the vorticity at heights z1 and z3 remains low for the vast majority of the radial domain. However, at z10 (bottom half of the conical section) it starts to rise dramatically at a normalized radius of r/R ≈ 0.5. This indicates that while the flow resembles the characteristics of an irrotational, free-vortex flow for the most part of the outer domain (r/R < 0.2) within the cylindrical and top part of the conical body, it becomes increasingly rotational at the bottom of the conical section, and specifically towards the apex. Furthermore, it is interesting to note that while the ratio of the increase of the vortex finder diameter for case B to case A is the same as that of case A to case C, i.e. (do)A / (do)B = (do)C / (do)A =2, the observed changes in the flow field characteristics are far greater in the former. In other words, the flow field undergoes more drastic changes when the vortex finder is doubled from case B to case A, than doubling the vortex finder diameter from case B to case C. This can imply that there is a critical value for do/du around which the vortex finder effect on the flow field is paramount, e.g. case B with do/du=1. However, when the do/du ratio exceeds that critical value, e.g. case C with do/du=4, the influence of changes in the vortex finder diameter on the flow field features dwindles. 114 Figure 5.3: Z-vorticity magnitude at z1=71.9 mm below the inlet, for vortex diameter cases. Figure 5.4: Z-vorticity magnitude at z3=110.1 mm below the inlet (top of the cone) for vortex diameter cases. 115 Figure 5.5: Z-vorticity magnitude at z10=377.3 mm below the inlet, for vortex diameter cases. The reason behind the remarkable dominance of the vorticity magnitude in case B at z1 over those of cases A and C is the fact that the underflow and overflow openings are of the same size, do/du=1. Due to the constant split ratio and mass influx across the cases, a much larger axial velocity field, and consequently the strongest rotational motion, i.e. vorticity, is generated in the case with the smallest do/du ratio (case B). In order to better explain the physics behind these observations the evolution of the vorticity magnitude across the hydrocyclone, for all cases A - C, is illustrated in Figure 5.6 -Figure 5.8). It can be seen that the overall patterns of the vorticity evolution in cases A and C are similar, as the flow is established inside the hydrocyclone. The vorticity profiles retain their shape and magnitude at heights inside the cylindrical body and top part of the conical body for these cases. However, further down the cone, and particularly in the vicinity of the apex, the rotationality rises considerably. For instance, at z10 the vorticity in case A is about 65%, and in case C about 35%, greater than that at z1, z4, and z8 in the same case. 116 Figure 5.6: Z-vorticity magnitude at various heights for case B. Figure 5.7: Z-vorticity magnitude at various heights for case A. 117 Figure 5.8: Z-vorticity magnitude at various heights for case C. The vorticity in case B, however, experiences a different evolution than that of cases A and C. It reaches its maximum inside the cylinder and in the proximity of the vortex finder. The vorticity declines as the fluid elements travel away from the vortex finder. Hence, the vorticity magnitude decreases despite the narrowing of the cross-sectional area in the cone. This indicates that the effect of the vortex finder opening size overrides that of the conservation of angular momentum due to the decrease in the cross-sectional area. In other words, the increase or preservation of the vorticity magnitude due to the cone angle is outgrown by the decrease in the vorticity magnitude due to the increased distance from the tip of the vortex finder. Thus, the net balance between these two effects results in a decline in the vorticity magnitude. It is observed that the magnitude drops about 40% from z1 to z8. At z10, however, the vorticity magnitude increases compared to z8. This is an indication that the proximity of the apex is far enough from the vortex finder for its effects to be diminished, and eventually be dominated, by the effect of the narrowing of the passage. The net balance between the two effects at this height leads to an increase in the vorticity, i.e. rotationality, compared to that at z8. 118 5.1.4 On Centrifugal Acceleration The centrifugal acceleration profiles at various heights are presented in Figure 5.9 - Figure 5.12). It is evident that the centrifugal acceleration in case B is greater, at all heights, than its counterparts in the other two cases. It is interesting to note that the centrifugal acceleration at z1 goes up by an order of magnitude from case C (do/du = 4) to case B (do/du = 1). Moreover, halving the ratio do/du in case B compared to case A leads to an eight-fold increase in the acceleration peak at z1. The same reduction in the ratio from case C to case A results in a two-fold increase in the peak. This can be an indication of the existence of a critical do/du ratio. Figure 5.9: Centrifugal acceleration vs. normalized radius at z1 for the vortex finder diameter cases. 119 Figure 5.10: Centrifugal acceleration vs. normalized radius at z4 for the vortex finder diameter cases. Figure 5.11: Centrifugal acceleration vs. normalized radius at z8 for the vortex finder diameter cases. The centrifugal acceleration increases, or preserves its magnitude, within the top half of the cone in cases A and C due to the reduction in the cross-sectional area and conservation of angular momentum. On the contrary, the acceleration decreases in case B at z4 and z8 compared to that at 120 z1. This is due to the dominance of the effect of far-away distance from the vortex finder over the effect of the reduction in the cross-sectional area. It is observed that at z8 and z9, the increase in the peak from case C to case A is more significant than that from case B to case A. This change in the relative discrepancy is attributed to the vanishing effect of the vortex finder diameter within the bottom half of the cone and towards apex. The acceleration in case B starts to increase at z9, which implies the area reduction becomes more influential than the vortex finder at this height (and further down towards the apex). Figure 5.12: Centrifugal acceleration vs. normalized radius at z9 for the vortex finder diameter cases. 5.1.5 On Bernoulli Gradient The Bernoulli gradient profiles for three different heights are presented in Figure 5.13 - Figure 5.15). It can be seen that at z1 the Bernoulli gradient in case B is an order of magnitude greater than that of cases A and C. This significant difference in the driving force of the flow field is responsible for the considerable discrepancy between the flow characteristics presented in the previous sections. The narrow opening of the vortex finder in case B leads to a dramatic increase in the Bernoulli function magnitude within a short amount of radial space, i.e. close proximity of the axis. This sharp increase, combined with the short radial distance, results in a far greater Bernoulli gradient in case B compared to the other cases. The distance from the vortex finder 121 weakens the Bernoulli gradient in case B. As a consequence, the discrepancy between peaks of the Bernoulli gradient profiles decreases at z8 (middle of the cone) compared to z1 (cylinder). Note that in cases A and C the effect of the narrowing of the passage is felt by the fluid elements earlier in the cone (at z8), and is manifested by the magnitude increase (Figure 5.14). The small vortex finder diameter in case B is suspected to delay such magnitude increase. As the fluid elements travel further down the cone and towards the underflow, e.g. at z9, the influence of the vortex finder vanishes and the reduction of the cross-sectional area dictates an increase in the Bernoulli gradient in case B. Figure 5.13: Bernoulli gradient vs. normalized radial position at z1 for cases A, B, and C. 122 Figure 5.14: Bernoulli gradient vs. normalized radial position at z8 for cases A, B, and C. Figure 5.15: Bernoulli gradient vs. normalized radial position at z9 for cases A, B, and C. 5.2 Effects of the Number of Inlets and Inlet Positioning The flow field features of two hydrocyclone geometries with dual inlets are presented in this subsection. A second inlet was added at a 180 circumferential degree with respect to the first inlet. 123 In one geometry (case D), the second inlet is at the same height as the first one, i.e. s = 0 with s being the difference between inlets axial positions. In the other geometry (case E), the second inlet is located below the first inlet by a distance equal to the inlet diameter, i.e. s = d. A schematic of these two dual-inlet geometries is shown in Figure 5.16. All other design and flow variables remain the same as the base case. The boundary condition and numerical methodology is the same as before. Note that the overall mass flow rate in these cases remains identical to that in the reference, single-inlet case. This means that the inlet velocity was adjusted such that the sum of mass influx through both inlets equal the mass flow rate in the single-inlet geometry. (a) (b) Figure 5.16: Geometries of hydrocyclone (a) case D with s=0 and (b) case E with s=d. 124 5.2.1 On Pressure The profiles of pressure drop, with respect to the inlet pressure, at various heights are illustrated in Figure 5.17. It is seen that the pressure drop behaves similarly for case D and the base, single-inlet case (case A). However, the pressure drop of case E is significantly lower, at all heights, than those of cases A and D. A second inlet at the same height as the first inlet appears to have little effect on the pressure drop. Figure 5.17: Pressure drop profiles at heights z1, z8, z9, z10; radial position normalized by the local wall radius. 125 Figure 5.18: Normalized pressure drop at heights z1, z8, z9, z10; radial position normalized by the local wall radius. Figure 5.18 shows the pressure drop profiles non-dimensionalized by the pressure difference between the overflow and underflow. It is interesting to note that the normalization causes the pressure drop profiles to coincide. 5.2.2 On Tangential Velocity The tangential velocity fields of cases A and D exhibit a similar trend, albeit discrepancies at their peak values. The tangential velocity in case E is evidently smaller than its counterparts. This is consistent with the observations for the pressure drop. A smaller pressure gradient, seen in Figure 5.16, generates a smaller tangential velocity field. 126 Figure 5.19: Tangential velocity profiles at heights, clockwise from top left, z1, z8, z9, z10; radial position normalized by the local wall radius. 5.2.3 On Centrifugal Acceleration Figure 5.20 compares the centrifugal accelerations for the dual-inlet cases with that of the base geometry. Again, cases A and D exhibit a similar pattern, with the acceleration peak slightly greater in case D. The centrifugal acceleration in case E is considerably lower than that in cases A and D. The lower pressure gradients in case D, and its smaller tangential velocity field observed in subsections 5.1.2 and 5.1.3, result in a lower centrifugal acceleration compared to cases A and E. 127 Figure 5.20: Centrifugal acceleration profiles at height z1 vs. normalized radial position. 5.3 Summary (i) It is observed that the tangential velocity is insensitive with changes in the vortex finder diameter size within the free-vortex region, i.e. for r/R > 0.5, where the flow is mainly irrotational. (ii) The tangential velocity peak is inversely proportional to the ratio of the size of the overflow (vortex finder) to the underflow openings, do/du. Decreasing do/du from 4 to 1 led to an approximately 100% increase in the tangential velocity peak. (iii) A non-dimensionlized pressure parameter was introduced for dual inlet cases which cause the pressure profiles of the single-inlet case and dual-inlet cases coincide irrespective of the inlet position. (iv) Lowering the axial position of the second inlet in the dual-inlet case resulted in a reduced centrifugal acceleration field, and consequently a decreased tangential velocity field. 128 Chapter 6: Sensitivity Study on the Reynolds Stress Model for the Hydrocyclone Flow Field 6.1 Axisymmetric Assumption The RANS equations coupled with the Reynolds stress turbulence model, i.e. equations (2.73) -(2.76), are numerically solved for an incompressible, unsteady, flow. The axisymmetric model has been chosen for this study due to its significantly lower computational cost compared to that of three-dimensional modeling. Generally speaking, the three-dimensional model produces more accurate results than does the axisymmetric model. However, it has been shown that the axisymmetric assumption gives reasonably accurate patterns for the hydrocyclone flow field (He et al.,1999; Hsieh, 1988). The flow asymmetry, i.e. three-dimensionality, is mainly caused by the hydrocyclone geometrical asymmetry due to the presence of a single inlet, narrowing passage in the apex vicinity and short-circuiting in the vortex finder proximity. Nonetheless, these asymmetry effects are confined to their respective regions: the proximities of the inlet, vortex finder and the apex. Velocity profiles in the main body, including the bottom of the cylindrical section and the top half of the cone, have been shown to be independent of the azimuthal angle (Hsieh, 1988). Therefore, the axisymmetry assumption is valid and justified to adopt for the purpose of this study. A three-dimensional investigation could be the subject of a future study, which may be performed using the insights into the coefficients’ influences, provided by the current study. 6.2 Computational Domain and Numerical Scheme Figure 6.1 shows the hydrocyclone grid and computational domain used in the simulations. Note that the depicted domain is mirrored around its axis, therefore it is twice as large as the simulated domain. The hydrocyclone dimensions are presented in Table 6.1. The simulations were performed by employing the Axisymmetric Swirl model in the CFD package ANSYS Fluent 16.01. The Semi-Implicit Method for Pressure Linked Equations (SIMPLE) was used for pressure-velocity coupling. The PREssure STagerring Option (PRESTO!), recommended for highly swirling flows (Ghadirian et al., 2013), was employed for pressure spatial discretization. The Quadratic Upwind Interpolation for Convective Kinematics (QUICK) algorithm was used for 129 the spatial discretization of the momentum, swirl and Reynolds stresses equations. This high order scheme was employed to minimize numerical diffusion and improve convergence (Ghadirian et al., 2013). Table 6.1: Hydrocyclone dimensions in mm. Labels have been illustrated in Figure 6.1. D din Lcyl Lcone Lv do du 102 19 102 405.7 62.5 12.7 6 Figure 6.1: The computational domain and mesh used in the simulations, mirrored with respect to the axis. The high concentration of mesh elements makes it difficult to clearly see the grid structure; thus, the zoomed-in cells are displayed for the sake of clarity. The geometrical axis of the hydrocyclone was specified as the axis of symmetry boundary condition. The pressure outlet boundary condition was applied at the overflow outlet with its value assigned/monitored to ensure a desired mass flow rate split. This boundary condition is justified as it is consistent with the experimental procedure of controlling the mass flow rate through the overflow; the split ratio in the LDV experiment was dictated by a valve opening. The underflow is assigned the pressure outlet boundary condition with the radial equilibrium pressure distribution, i.e. radial gradient allowed, which is a more realistic condition considering the narrowing of the passage towards the apex. The impermeable, no-slip wall boundary condition is used for the hydrocyclone roof, the vortex finder tube, and the cylindrical and conical walls. The velocity 130 boundary condition is considered for the inflow. The flow enters a conventional hydrocyclone through a single tube at one side of the cyclone axis. However, an axisymmetric model assumes the inlet is continuous along the circumference of the cyclone, i.e. starting from θ=0 at one end and connecting at θ=360 at the other end. In other words, the flow enters the cyclone through a complete circumferential ring, referred to as the ring inlet. The ring inlet in axisymmetric models is an equivalent of a single tube in three-dimensional models and demands a modified full 360° inlet boundary condition. It is imperative to ensure that the resulting flow rate is the same as that in the experiment, both in terms of mass inflow and momentum influx. Therefore, the radial velocity component is responsible for convecting mass into the hydrocyclone, while the tangential velocity component produces the equivalent swirl at the inlet. The following relationships have been suggested in the literature for the inlet radial and swirl velocity components (Malhotra et al., 1994; He et al., 1999; Hsieh, 1988), 𝑈𝑟 =𝑄𝑖𝑛𝜋𝐷𝑑𝑖𝑛 (6.1) 𝑈𝜃 =4𝑄𝑖𝑛𝜋𝑑𝑖𝑛2 (6.2) where 𝑄𝑖𝑛, 𝐷, and 𝑑𝑖𝑛 are the inlet volumetric flow rate, hydrocyclone diameter, and the inlet tube diameter, respectively. The principle behind the ring inlet assumption is that the flow entering through a single-tube entry in conventional cyclones will quickly form a swirling layer around the inside of the cylinder wall. Thus, the inlet velocity is transformed into swirl at the intersection of the inlet and the cylindrical section. This result mimics that of the ring inlet in the axisymmetric model. The velocity components at the inlet were defined to result in a mass influx of 0.72 kg/s, corresponding to the experimentally measured velocity at the inlet. The hydraulic diameter was defined the same as the inlet diameter with a turbulent intensity of 10%. The k-ε model was used for initiation, after which the LPS was chosen to advance several more time steps. Once a crude flow field was established, i.e., with the pressure and momentum in balance, the turbulence model was switched to the QPS. This initialization process was to avoid numerical instability and convergence difficulties caused 131 by the QPS. An overall time step size of 5×10-4 s was used. The computational time was approximately 168-240 hours, for each case, using 12 cores on an i7 CPU machine. 6.3 Grid Sensitivity The mesh independency study was performed using the benchmark turbulence model with the default values for the RSM constant, and a mass influx of 0.72 kg/s (velocity of 2.60 m/s at the inlet) with a split ratio of 50%. Figure 6.2 shows the swirl velocity extracted at an axial location within the cone for three different mesh sizes; coarse, intermediate, and fine grids with approximately 9000, 20,000 and 36,000 mesh elements. The velocity field comparison shows a maximum of 2% discrepancy between the predictions from the intermediate mesh and fine mesh results. A similar trend exists across the hydrocyclone. Therefore, the rest of the simulations were carried out using the fine grid. Figure 6.2: Swirl velocity profile for different grid sizes at z = 364.5 mm. 6.4 Experimental Validation The numerical results were compared with LDV measurements (Mackenzie 2014) for a hydrocyclone with identical geometrical and flow parameters. It is worth mentioning that the error associated with the LDV measurements was reasonably small, with a maximum error of 10% for 132 the tangential velocity field at a location near the axis. For more details on the experimental measurements, including the error analysis as well as the flow loop, see Mackenzie (2014) and MacKenzie et al. (2014). Figure 6.3 compares tangential velocity profiles calculated from the CFD simulations for both the LPS and QPS models, against the LDV data. The elevation presented in Figure 6.3(a) is in the cylindrical body of the cyclone, whereas radial positions in Figure 6.3(b) - (d) are located in the conical section. It can be observed that while the LPS produces a closer match to the experimental data within the free-vortex region, it considerably under-predicts the swirl velocity within the forced-vortex core. However, the QPS results produce closer predictions to the LDV data within the rigid body rotation core, despite some discrepancy between the peak values. Thus, in spite of some improvements within the forced-vortex core, the QPS has not fully resolved the under-prediction problem, as is evident in Figure 6.3(a) and (d). Moreover, the QPS led to slight over-predictions within the free-vortex region. Figure 6.3: Swirl velocity profiles at (a) z=97.1 mm, (b) z= 115.4 mm, (c) z=135.5 mm, and (d) z=154.5 mm, away from the cyclone roof; radial position is normalized by the local radius. 133 The axial velocity profiles extracted at different axial locations are plotted against the LDV data in Figure 6.4. It is seen that the predictions of the two CFD models for the axial velocity profiles within the outer, free-vortex region, which approximately amounts to 80% of the radial domain, of the hydrocyclone closely agree with the experimental measurements. Figure 6.4: Axial velocity profiles at (a) z=97.1 mm, (b) z= 115.4 mm, (c) z=135.5 mm, and (d) z=154.5 mm, away from the hydrocyclone roof; radial position is normalized by the local radius. However, within the inner core, i.e. r/R < 0.2, the predictions of both the LPS and QPS start to deviate from the experimental data. It is evident from Figure 6.4 that both models over-predict the peak of the axial velocity at all heights. The discrepancy grows larger as the cyclone axis is approached; the LPS model leads to the most substantial deviation from the experimental data. Therefore, the QPS model gives a better performance than the LPS model. However, the over-prediction problem persists in the QPS model as well, even though to a lesser extent than the LPS. 134 Figure 6.5: Swirl velocity profiles at (a) z=211.8 mm, (b) z= 250 mm, (c) z=370.3 mm, and (d) z=364.5 mm, away from the cyclone roof; radial position is normalized by the local radius. 6.5 Effect of RSM Coefficients The twelve coefficients of the quadratic pressure-strain RSM were increased, one at a time, to double their default numerical values in the default case (case 0). This results in 12 different cases, called doubled cases hereafter, (cases 1, 3, 5 … 23), in order to investigate the effect of increasing every single constant on the model’s performance. A similar process was completed by reducing the coefficients to half their default values, referred to as halved cases (cases 2, 4, 6 … 24), in order to study the influence of decreasing the constants. This results in 24 different cases in addition to the default case (case 0). The response of the model to these changes is presented in the following sections. The cases are categorized, based on their response to their coefficient value change, into three different groups, insensitive, moderately sensitive, and highly sensitive. 135 6.5.1 Doubled Cases Table 6.2 presents the coefficients values used in doubled cases; in each case only the value of one constant is increased at a time. For example, in case 1 Cμ is doubled from its default value of 0.09 to 0.18 while the rest of the coefficients maintain their default values, i.e. identical to those of the default case 0. The changed values for each case are highlighted in Table 6.2. Note that cases 3, 11, 13, and 21 were the only cases for which their respective constants were increased to less than double the default values. This was owing to convergence difficulties encountered when doubling their values. The swirl velocity profiles extracted at an axial location inside the cylindrical body and an axial location within the conical section are illustrated in Figure 6.6 and Figure 6.7). The axial locations are stated with reference to their distance from the cyclone roof (top wall). The response to the changes in the RSM coefficients values was observed to follow a similar trend at other axial heights. Table 6.2: RSM coefficients used in doubled cases; those highlighted in green are double the default values in case 0. Cmu C1-Epsilon C2-Epsilon C1-SSG-PS C1'-SSG-PS C2-SSG-PS C3-SSG-PS C3'-SSG-PS C4-SSG-PS C5-SSG-PSTKE Prandtl TDR Prandtl Case 0 0.09 1.44 1.83 3.4 1.8 4.2 0.8 1.3 1.25 0.4 1 1.3Case 1 0.18 1.44 1.83 3.4 1.8 4.2 0.8 1.3 1.25 0.4 1 1.3Case 3 0.09 1.54 1.83 3.4 1.8 4.2 0.8 1.3 1.25 0.4 1 1.3Case 5 0.09 1.44 3.66 3.4 1.8 4.2 0.8 1.3 1.25 0.4 1 1.3Case 7 0.09 1.44 1.83 6.8 1.8 4.2 0.8 1.3 1.25 0.4 1 1.3Case 9 0.09 1.44 1.83 3.4 3.6 4.2 0.8 1.3 1.25 0.4 1 1.3Case 11 0.09 1.44 1.83 3.4 1.8 8 0.8 1.3 1.25 0.4 1 1.3Case 13 0.09 1.44 1.83 3.4 1.8 4.2 1 1.3 1.25 0.4 1 1.3Case 15 0.09 1.44 1.83 3.4 1.8 4.2 0.8 2.6 1.25 0.4 1 1.3Case 17 0.09 1.44 1.83 3.4 1.8 4.2 0.8 1.3 2.5 0.4 1 1.3Case 19 0.09 1.44 1.83 3.4 1.8 4.2 0.8 1.3 1.25 0.8 1 1.3Case 21 0.09 1.44 1.83 3.4 1.8 4.2 0.8 1.3 1.25 0.4 2 1.3Case 23 0.09 1.44 1.83 3.4 1.8 4.2 0.8 0.65 1.25 0.4 1 2136 Figure 6.6: Swirl velocity profiles at z=85.4 mm below the roof (in the cylinder) for the QPS cases with increased constants vs. the default case. Figure 6.7: Swirl velocity profiles at z=225.40 mm below the roof (in the cone) for the QPS cases with increased constants vs. the default case. 137 6.5.1.1 Insensitive The swirl velocity profiles of cases 11 and 17 coincide to those of case 0. Therefore, increasing their respective constants caused little to no changes in the swirl velocity, and consequently in the centrifugal field. The same behaviour was observed at other axial locations. Thus, the C2-SSG-PS and C4-SSG-PS coefficients are classified as insensitive in regard to a value increase as high as 100%. 6.5.1.2 Moderately Sensitive The centrifugal fields predicted by cases 1, 7, 9, and 19 are very similar to that of the default case. It is evident from Figure 6.6 and Figure 6.7) that the swirl velocity profiles of these cases closely match that of case 0 for the vast majority of the domain. However, there is a slight difference between the predictions in the vicinity of the intersection of the free-vortex and forced-vortex regions, including the maximum swirl velocity. In all four cases, the increase of their respective coefficients resulted in a slight decline of the swirl velocity at its peak, as well as within a distance of approximately 5% of the radial position corresponding to the maximum swirl on both sides of the peak. It is worth mentioning that this discrepancy appears to grow larger as the flow approaches the apex. Considering the relatively minimal changes produced in these cases, their respective coefficients, Cμ, C1-SSG-PS, C1′-SSG-PS, and C5-SSG-PS are categorized as moderately sensitive with respect to increasing their values (up to a 100% increase). 6.5.1.3 Highly Sensitive The deviation from the default tangential velocity field is rather significant in the results produced by cases 3, 5, 13, 15, 21, and 23. The swirl velocity profiles obtained from cases 5 and 21 exhibit the most dramatic declines compared to case 0; they fail to even capture the distinct pattern of the Rankine vortex, the main characteristic of the hydrocyclone flow field. The swirl velocities obtained by these two cases are somewhat flat with almost no features that could suggest a combination of a free and forced vortex. This is an adverse change since it worsens the under-prediction issue associated with the RSM model with its default coefficients. On the contrary, the results of case 15 show an improvement both in terms of the maximum velocity and the better 138 agreement for almost half of the free-vortex region. However, this improvement is impaired by the considerable fall (as high as 50%) of the swirl velocity profiles in the majority of the forced-vortex core as well as the inner half of the free-vortex region. Cases 3, 13, and 23 produced a sharp increase in the peak swirl velocity, ranging from 42% to 57% increase in the maximum swirl velocity predictions relative to the default case. Moreover, note that the changes in the free-vortex region are minimal; the predictions of cases 3, 13, and 23 match that of the default case in approximately 75% to 90% of the outer free-vortex core. This dramatic increase in the peak velocity and within the rigid-body rotation core, along with the relatively unchanged free-vortex region, is promising. It suggests that fine-tuning the respective constants could lead to a set of new model constants which helps overcome the under-prediction problem of the RSM. Based on these observations coefficients Cε1, Cε2, C3-SSG-PS, C3′-SSG-PS, TKE Prandtl Number and TDR Prandtl Number are categorized as highly sensitive with respect to the increases in their values. It should be noted that Cε1, C3-SSG-PS, and TDR Prandtl Number resulted in desirable responses and therefore should be investigated further. 6.5.2 Halved Cases The RSM coefficients used for the halved cases are presented in Table 6.3. Similarly one coefficient at a time is changed to half its default value in case 0, while the remaining eleven coefficients retain their default values. A similar analysis is performed for the halved cases, which is presented below. The swirl velocity profiles are presented in Figure 6.8 andFigure 6.9); the axial locations are the same as those used for the doubled cases. 139 Table 6.3: RSM coefficients used in halved cases; those in green are half the default values in case 0. Figure 6.8: Swirl velocity profiles at z=85.4 mm below the roof for the QPS cases with reduced constants vs. the default case. Cmu C1-Epsilon C2-Epsilon C1-SSG-PS C1'-SSG-PS C2-SSG-PS C3-SSG-PS C3'-SSG-PS C4-SSG-PS C5-SSG-PSTKE Prandtl TDR Prandtl Default 0.09 1.44 1.83 3.4 1.8 4.2 0.8 1.3 1.25 0.4 1 1.3Case 2 0.045 1.44 1.83 3.4 1.8 4.2 0.8 1.3 1.25 0.4 1 1.3Case 4 0.09 0.77 1.83 3.4 1.8 4.2 0.8 1.3 1.25 0.4 1 1.3Case 6 0.09 1.44 0.915 3.4 1.8 4.2 0.8 1.3 1.25 0.4 1 1.3Case 8 0.09 1.44 1.83 2.45 1.8 4.2 0.8 1.3 1.25 0.4 1 1.3Case 10 0.09 1.44 1.83 3.4 0.9 4.2 0.8 1.3 1.25 0.4 1 1.3Case 12 0.09 1.44 1.83 3.4 1.8 2.1 0.8 1.3 1.25 0.4 1 1.3Case 14 0.09 1.44 1.83 3.4 1.8 4.2 0.4 1.3 1.25 0.4 1 1.3Case 16 0.09 1.44 1.83 3.4 1.8 4.2 0.8 0.65 1.25 0.4 1 1.3Case 18 0.09 1.44 1.83 3.4 1.8 4.2 0.8 1.3 0.625 0.4 1 1.3Case 20 0.09 1.44 1.83 3.4 1.8 4.2 0.8 1.3 1.25 0.2 1 1.3Case 22 0.09 1.44 1.83 3.4 1.8 4.2 0.8 1.3 1.25 0.4 0.5 1.3Case 24 0.09 1.44 1.83 3.4 1.8 4.2 0.8 0.65 1.25 0.4 1 0.650140 Figure 6.9: Swirl velocity profiles at z=225.4 mm below the roof for QPS cases with reduced constants vs. the default case. 6.5.2.1 Insensitive Cases 2, 8, and 12, appear to show almost zero sensitivity to decreases in the numerical values of their respective coefficients. Their swirl velocity profiles are observed to collapse onto that of case 0. The same insensitivity exists at other axial locations including near the apex. Although the swirl profiles of cases 18 and 20 do not perfectly match the default case results, particularly in a small proximity of the peak, the difference is insignificant. These observations led us to categorize the corresponding coefficients, i.e. Cμ, C1-SSG-PS, C2-SSG-PS, C4-SSG-PS and C5-SSG-PS, as insensitive to a decrease of 100% in their numerical values. 6.5.2.2 Moderately Sensitive The swirl velocity profiles obtained from cases 6 and 10 closely follow the swirl velocity of case 0 all along the radial direction at an axial location of 85.4 mm. However, their peak values go up slightly compared to the maximum swirl of case 0. Thus, the coefficients Cμ and C1′-SSG-PS are categorized as moderately sensitive to a decrease of 100% in their numerical values. Note that 141 these changes are desirable responses which can eliminate (or at least alleviate) the RSM under-prediction problem. 6.5.2.3 Highly Sensitive Cases 4 and 14 responded dramatically to the reduction of their coefficients to half their original numerical values. They predicted an almost constant swirl velocity profile with no distinctive free or forced vortex regions. These changes, though substantial, are undesirable and only deteriorate the RSM under-prediction problem. The swirl velocity of case 24 closely matches its default counterpart throughout the free-vortex region. However, there are sharp declines in the forced-vortex core and in the peak magnitude. The peak appears to be broader as well. Considerable but desirable changes are observed in the swirl velocity predictions of cases 16 and 22. It has been illustrated in Figure 6.8 that their profiles closely match that of the default case all along the free-vortex region. However, there is a significant increase in their predictions within the forced-vortex region and in particular for the peak swirl magnitude. This increase is a desirable response which should be investigated more closely. Thus, coefficients Cε1, C3-SSG-PS, C3′-SSG-PS, TKE Prandtl Number and TDR Prandtl Number are categorized as highly sensitive to a decrease of 100% in their numerical values. It should be noted that only C3′-SSG-PS and TKE Prandtl Number exhibited a positive and desirable response to the reduction in their numerical value. 6.6 Hydrodynamics of the Flow Field The radial velocity component plays a crucial role in hydrocyclone particle separation. In this section, we examine the radial momentum equation more closely. To do so, the terms in the radial momentum equation, equation (2.74), are computed separately to quantify their contribution to the force balance at various spatial locations in the hydrocyclone. The flow variables, including the velocity components, pressure and the relevant Reynolds stress components, were extracted and the derivatives were numerically calculated based on cell values. The balance between the different 142 forces is illustrated in Figure 6.10. Note that the left-hand side of equation (2.74) is referred to as the convective forces and the pressure gradient term on the right-hand side is denoted by pressure force. The second term on the right-hand side of equation (2.74) represent the viscous forces, and turbulence force denotes the sum of the Reynolds stresses terms (the last three terms on the right-hand side of the radial momentum equation). It is evident from Figure 6.10 (left) that the convective and pressure forces are by far the most dominant forces in the radial momentum equation. For the sake of illustration, the viscous and turbulence forces have been depicted on a separate graph, Figure 6.10 (right). It can be seen from this graph that the viscous forces are several orders of magnitude smaller than the convective and pressure forces throughout the domain. Similarly, turbulence forces are at least three orders of magnitudes smaller than the pressure and convective forces. However, one should note that despite the significantly smaller contribution of the Reynolds stresses to the force balance, they cannot and should not be neglected in the momentum equation due to their significance in the turbulence generation and sustenance. Their presence is imperative in establishing the oscillatory and highly swirling flow field within the hydrocyclone. Figure 6.10: Profiles of all forces in equation (2.74) (left), and profiles of viscous and turbulence forces (right) at z=102 mm from the cyclone roof. These results correspond to the QPS model. 143 6.7 Force Balance Analysis for LPS and QPS The convective and pressure forces for the LPS and QPS models are presented in Figure 6.11. It is observed that in both models the convective and pressure forces match one another almost perfectly. Therefore, they essentially govern the momentum equation in the radial direction and particularly so within the forced-vortex core. Interestingly, while these two forces have a similar effect and order of magnitude in the LPS and QPS models, their patterns appear to be drastically different. While the pressure and convective forces predicted by these two models match one another closely in the free-vortex zone (specifically in areas with r/R > 0.2), they deviate significantly from one another in the rigid-body-rotation core (r/R < 0.1). The pressure force and convective forces in the LPS model suggest a similar trend to that of the swirl velocity; starting from low values in the outer, free-vortex region, reaching its peak and then plummeting in the inner core. This is reminiscent of the Rankine vortex pattern of the swirl velocity. The QPS model, however, predicts an opposite trend in the forced-vortex region. The QPS model results indicate that the pressure and convective forces continue growing larger in the vicinity of the centre with no sign of any decline in their magnitude, unlike the LPS model. We suspect this stark difference in the trend and magnitude of the dominant forces in the LPS and QPS can be responsible for the discrepancy in the LPS and QPS predictions observed both in our simulations and elsewhere (Cullivan et al., 2003a; Cullivan et al., 2004). 144 Figure 6.11: Pressure and convective forces profiles for the QPS and LPS models at z=102 mm away from the roof - the border of the conical and cylindrical bodies. Figure 6.12: The ratio of the centrifugal force to the convective forces in the QPS model, at z=102 mm from the roof. Figure 6.12 shows the ratio of the centrifugal force term, −𝜌𝑈𝜃2𝑟, to the convective forces, i.e. the sum of the three terms on the left hand side of equation (2.74), for the QPS model at the border of 145 the cylindrical and conical bodies, z=102 mm. It is observed that the centrifugal force is by far the dominant term amounting to almost 100% of the convective forces throughout the domain. Its contribution drops only by a maximum of 5%, to 95%, in the vicinity of the wall. This highlights the overwhelming significance of the centrifugal force compared to the other two convective forces and suggests that they can be neglected without significant error. Note that this result can support the assumptions made in analytical literature where the r-momentum was simply reduced to the balance between the pressure gradient and the centrifugal force, i.e. 𝜌𝑈𝜃2𝑟=𝜕𝑃𝜕𝑟 (Rietema and Maatschappij, 1961). 6.8 Summary CFD simulations of a 102 mm diameter hydrocyclone were performed using the commercial package ANSYS Fluent 16.01. An axisymmetric model was employed along with the Reynolds stress turbulence model. The CFD predictions were in good agreement with experimental LDV data. A numerical study was conducted to investigate the influence of the default constants of the standard RSM on the model’s predictions. One constant was changed at a time while the rest retained their standard values. Furthermore, a force balance analysis was performed and the contributions of different forces in establishing the flow field pattern were quantified. The possible hydrodynamic reason behind the observed discrepancy in the predictions by the LPS and QPS, often reported in literature, was explained. The following conclusions can be drawn: The predictions from the QPS model were a closer match to the experimental data than those of the LPS model. However, both models under-predicted the swirl velocity field, particularly the peak where the free and force vortices meet, at all heights. The dominance of convective and pressure forces in the radial momentum equation was quantified. They are shown to be at least three orders of magnitudes greater than the turbulence force and six to seven orders of magnitudes greater than the viscous forces. Hence, the radial momentum equation could well be treated as a balance between the pressure and convective forces. The sum of all convective forces is almost equal to the centrifugal force, 𝜌𝑈𝜃2𝑟, for approximately 80% of the domain. In the remaining 20% of the domain, including the wall 146 proximity, the centrifugal force contributes as high as 95% to the total amount of the convective forces. Therefore, by reducing the left-hand side of the radial momentum equation (2.73), sum of the convective forces, to the centrifugal force only, 𝜌𝑈𝜃2𝑟, one would at worst introduce a maximum error of 5%. Turbulence is observed to reach its peak within the forced-vortex core; for the QPS model at the cyclone axis and for the LPS model slightly away from the axis, and before the position of the maximum swirl. There was little difference between the behaviour of the pressure and convective forces in the LPS and QPS for approximately 90% of the domain, which includes the entire outer free-vortex region. In the inner forced-vortex core, the behaviour of the pressure and convective forces in the LPS deviates radically from that in the QPS. The pressure and convective forces increase dramatically in the inner core in the QPS, whereas their counterparts in the LPS decline sharply. The drastically different trends can be responsible for different predictions by the LPS and QPS. Changes in the values of Cε1, Cε2, C3-SSG-PS, C3′-SSG-PS, TKE Prandtl, and TDR Prandtl coefficients prompted the greatest, but not necessarily desirable, changes in CFD predictions of the swirl velocity. Desirable responses were elicited from the model by increasing Cε1, C3-SSG-PS, and TDR Prandtl as well as decreasing Cε2, C3′-SSG-PS, and TKE Prandtl. These changes all resulted in greater swirl velocity peaks. This can be insightful in fine-tuning the model constants to resolve the RSM under-prediction issue. Changes in the respective values of C2-SSG-PS and C4-SSG-PS produced no change in the CFD predictions compared to those of the benchmark, default case. Similarly, the simulation results exhibit little-to-no changes with decreasing the values of Cμ, C1-SSG-PS, C5-SSG-PS. 147 Chapter 7: Summary and Conclusions 7.1 Objective 1 Our first objective was to gain a deeper understanding of fluid forces in swirling flows. We found a critical threshold that determines when convective forces dominate over turbulence forces. This threshold determination is based on an “integral measure criterion” of local forces in the radial and axial directions. The threshold itself is defined by a dimensionless number, N, based on the ratio of the circumferential and axial flow velocities. The critical value was found to be Ncr =0.45. Above this, convection dominates; below it, turbulence dominates. This finding will facilitate selection of CFD models in future work to optimize cost and accuracy for modelling swirling flows. For example, k-ε models suffice when Ncr < 0.45, but more complex models are required for higher values. 7.2 Objective 2 The aim of this work was to study, in detail, the effect of cone angle on the hydrocyclone flow hydrodynamics. We found a new dimensionless pressure parameter to determine local pressures along the axis of hydrocyclones. This dimensionless number is the ratio of two Euler numbers based on the inflow, overflow, underflow, and local pressures and velocities. It was shown that the dimensionless pressures, for three cone angles tested, all fell on a common curve. This finding will assist in scaling hydrocyclones with differing cone angles. 7.3 Objective 3 Our third objective was to investigate the default constants of the standard quadratic pressure-strain (QPS) Reynolds Stress Model. We changed one constant at a time, and from the responses we classified the constants into strong, medium, and zero influence. In doing so, we identified candidate constants for fine tuning RSM as a hydrocyclone-specific turbulence model. Desirable responses were elicited from the model by increasing Cε1, C3-SSG-PS, and TDR Prandtl as well 148 as decreasing Cε2, C3′-SSG-PS, and TKE Prandtl. These changes all resulted in greater swirl velocity peaks. These findings will assist in fine-tuning model constants to resolve the RSM under-production issue. Furthermore, based on these findings, we speculate that the hydrodynamic reasons for previously reported discrepancies between predictions from the QPS and LPS models is that pressure and convective forces increase at the inner core of the QPS model, whereas they decrease sharply in the LPS model. 149 Chapter 8: Implications and Limitations of the Research and Recommendations for Future Work 8.1 Objective 1 Future work could build upon these findings to develop means to assess cost-effective yet accurate CFD models for other geometries. As an example, the flow field at low rotation numbers (N < Ncr) can be simulated and compared to experimental data. A further extension could be to extend and test the critical number to encompass three-dimensional models. The large simulations matrix led us to use axisymmetric model for this part of the study. Despite the reasonable accuracy shown for the axisymmetric predictions, a three-dimensional approach would have been the most accurate model to use. 8.2 Objective 2 The proposed dimensionless pressure group should be tested for hydrocyclones with different geometrical features to test its universality. Furthermore, the dimensionless number can be tested for the same geometry with different operational parameters such as throughput and split ratio. Similar to objective 1, a limitation of this part of the study was the use of axisymmetric modelling. Despite reasonably accurate predictions of the model, validated with the experimental data, a 3D modelling would provide the most accurate predictions for the hydrocyclone flow field. 8.3 Objective 3 Building upon the results of this work, future work can aim at the ultimate goal of a hydrocyclone-specific turbulence model. A useful next step would be the use of the Design of Experiment (DOE) approach instead of the “One-variable-at-a-time” approach employed in this work. Success with this would greatly reduce the scope of future experiments necessary to verify CFD models. 150 In addition, a similar OVAT technique (or DOE) should be used for a number of hydrocyclone geometries to investigate possible geometry-dependence of the outcome performance to changes in the input variables, i.e. the RSM constants. This will be a necessary step in developing a geometry-independent turbulence model for hydrocyclones. Eventually, the same approach (OVAT or a newly planned DOE) should be applied to three-dimensional hydrocyclone models to investigate the extent of, and eliminate, possible effects of the axisymmetric assumption on the outcome performance. This axisymmetric study serves as an intermediate step taken towards the ultimate goal of a 3D geometry-independent, hydrocyclone-specific turbulence model. The simulated cases were without an air core. This was due to the fact that no air core was observed in experiments for an identical hydrocyclone with the same conditions as those used in the simulations. However, multiphase modelling, VOF (Volume of Fluid) or Mixture Model, could be employed to resolve an air core. Another limitation was a lack of particle separation study. One-way coupling of DPM (Discrete Phase Model) or two-way coupling, should be used to investigate particle motion and the separation mechanism. It should be noted the existence of a high solid concentration in the flow, such as pulp slurry with medium to high pulp consistency, can lead to rheological change in the flow. Turbulence will be dampened, the tangential velocity diminishes, and the free-vortex region starts to vanish at high solid concentrations. This was another limitation of the work as the presented results in this thesis may be applicable only to dilute slurries and suspensions where the solid particles do not disturb the fluid flow significantly; particles Stokes number is low. While being a limitation, this effect can be a future work goal as well. 151 References Ansys, Inc. 2011. “ANSYS Fluent Theory Guide.” Canonsburg, Pennsylvania 15317 (November): 794. doi:10.1016/0140-3664(87)90311-2. Averous, J., and Fuentes, R., 1997. “Advances in the Numerical Simulation of Hydrocyclone Classification.” Canadian Metallurgical Quarterly 36 (5): 309–14. doi:10.1179/cmq.1997.36.5.309. Berges, A. 1935. British Patent No. 455, 845, issued 1935. Bergström, J., 2006. “Flow Field and Fibre Fractionation Studies in Hydrocyclones”, PhD thesis, Royl Institute of Technology (KTH), Stockholm, Sweden. https://www.diva-portal.org/smash/get/diva2:11075/FULLTEXT01.pdf. Bergstrom, J., Vomhoff, H. and Soderberg, D., 2007. “Tangential Velocity Measurements in a Conical Hydrocyclone Operated with a Fibre Suspension.” Minerals Engineering 20 (4): 407–13. doi:10.1016/j.mineng.2006.09.008. Boysan, F., Ayers, W. H. and Swithenbank, J., 1982. “A FUNDAMENTAL MATHEMATICAL-MODELING APPROACH TO CYCLONE DESIGN.” Transactions of The Institution of Chemical Engineers 60 (4): 222–30. Bradley, D., Pulling, D. J. 1959. “Flow Pattern in the Hydraulic Cyclone and Their Interpretation in Terms of Performance.” Trans Inst Chem Eng 37: 34–45. Bradley, D. 1965. The Hydrocyclone. 1st ed. Pergamon Press Ltd. Brennan, Mathew S. 2006. “CFD Simulations of Hydrocyclones with an Air Core: Comparison Between Large Eddy Simulations and a Second Moment Closure.” Chemical Engineering Research and Design 84 (6): 495–505. doi:10.1205/cherd.05111. Brennan, M., Holtham, P., and Narasimha, M., 2009. “CFD Modelling of Cyclone Separators : Validation Against Plant Hydrodynamic Performance.” In Seventh International 152 Conference on CFD in the Minerals and Process Industries. CSIRO, Melbourne, Australia. Bretney, E. 1891. The Hydrocyclone. 453, issued 1891. Chen, C. J., and Jaw, S. Y., 1998. Fundamentals of Turbulence Modeling. Washington, DC: Taylor & Francis. Chen, W., Zydek, N., and Parma, F., 2000. “Evaluation of Hydrocyclone Models for Practical Applications.” Chemical Engineering Journal 80 (1–3): 295–303. doi:10.1016/S1383-5866(00)00105-2. Chiné, B., and Concha, F., 2000. “Flow Patterns in Conical and Cylindrical Hydrocyclones.” Chemical Engineering Journal 80 (1–3): 267–73. doi:10.1016/S1383-5866(00)00101-5. Chu, L. Y., Chen, W. M., and Lee, X. Z., 2000. “Effect of Structural Modification on Hydrocyclone Performance.” Separation and Purification Technology 21 (1–2): 71–86. doi:10.1016/S1383-5866(00)00192-1. Concha, F. 2007. “Flow Pattern in Hydrocyclones.” KONA 25: 97–132. Cullivan, J. C., Williams, R. A. and Cross, C. R., 2003a. “Understanding the Hydrocyclone Separator through Computational Fluid Dynamics.” Trans IChemE 81 (April): 455–66. doi:10.1205/026387603765173718. ———. 2003b. “New Insights into Hydrocyclone Operation.” Particulate Science and Technology 21 (1). Informa UK Ltd: 83–103. doi:10.1080/02726350307500. Cullivan, J. C., Williams, R. A., Dyakowski, T., and Cross, C. R., 2004. “New Understanding of a Hydrocyclone Flow Field and Separation Mechanism from Computational Fluid Dynamics.” Minerals Engineering 17 (5): 651–60. doi:10.1016/j.mineng.2004.04.009. Dabir, B., Petty, C. A., 1986. “Measurement of Mean Velocity Profiles Using Laser Doppler Anemometry.” Chemical Engineering Communications 48 (4–6): 377–88. Dai, G. Q., Li, J. M., and Chen, W. M., 1999. “Numerical Prediction of the Liquid Flow within a 153 Hydrocyclone.” Chemical Engineering Journal 74 (3): 217–23. doi:10.1016/S1385-8947(99)00044-3. Daly, B. J., and Harlow, F. H., 1970. “Transport Equations in Turbulence.” Physics of Fluids 13 (11). American Institute of PhysicsAIP: 2634. doi:10.1063/1.1692845. Deardorff, J. W. 1970. “A Numerical Study of Three-Dimensional Turbulent Channel Flow at Large Reyn,olds Numbers.” Journal of Fluid Mechanics 41 (2). Cambridge University Press: 453. doi:10.1017/S0022112070000691. Delgadillo, J. A., and Rajamani, R. K., 2009. “Computational Fluid Dynamics Prediction of the Air-Core in Hydrocyclones.” International Journal of Computational Fluid Dynamics 23 (2): 189–97. doi:10.1080/10618560902724893. Delgadillo, J. A., and Rajamani, R. K., 2005. “A Comparative Study of Three Turbulence-Closure Models for the Hydrocyclone Problem.” International Journal of Mineral Processing 77 (4): 217–30. doi:10.1016/j.minpro.2005.06.007. ———. 2007. “Exploration of Hydrocyclone Designs Using Computational Fluid Dynamics.” International Journal of Mineral Processing 84 (1–4): 252–61. doi:10.1016/j.minpro.2006.07.014. Dyakowski, T., and Williams, R. A., 1993. “Modelling Turbulent Flow within a Small-Diameter Hydrocyclone.” Chemical Engineering Science 48 (6): 1143–52. doi:10.1016/0009-2509(93)81042-T. Feiz, A A, Ould-Rouis, M., and Lauriat, G, 2003. “Large Eddy Simulation of Turbulent Flow in a Rotating Pipe.” International Journal of Heat and Fluid Flow 24 (3): 412–20. doi:10.1016/S0142-727X(02)00241-2. Ferziger, J. H. 1996. “Large Eddy Simulation.” In Simulation and Modeling of Turbulence Flows, edited by B. Thomas Gatski, Yousuff M. Hussaini, and John L. Lumley. New York. Fisher, M. J., and Flack, R. D., 2002. “Velocity Distributions in a Hydrocyclone Separator.” 154 Experiments in Fluids 32 (3): 302–12. doi:10.1007/s003480100344. Freeman, H. 1936. U. S. Patent No. 2,102,525, issued 1936. Fu, S., Launder, B. E., and Leschziner, M. A., 1987. “Modeling Strongly Swirling Recirculating Jet Flows with Reynolds Stress Transport Closure.” In Sixth Siyposium on Turbulent Shear Flows. Toulouse, Fracne. Germano, M., Piomelli, U., Moin, P., and Cabot, W. H., 1991. “A Dynamic Subgrid‐scale Eddy Viscosity Model.” Physics of Fluids A: Fluid Dynamics 3 (7). American Institute of PhysicsAIP: 1760–65. doi:10.1063/1.857955. Ghadirian, M., Hayes, R. E., Mmbaga, J., Afacan, A., and Xu, Z., 2013. “On the Simulation of Hydrocyclones Using CFD.” Canadian Journal of Chemical Engineering 91 (5): 950–58. doi:10.1002/cjce.21705. Ghodrat, M., Kuang, S. B., Yu, A. B., Vince, A., Barnett, G. D., and Barnett, P. J., 2014a. “Numerical Analysis of Hydrocyclones with Different Vortex Finder Configurations.” Minerals Engineering 63: 125–38. doi:10.1016/j.mineng.2014.02.003. Ghodrat, M., Kuang, S. B., Yu, A. B., Vince, A., Barnett, G. D., and Barnett, P. J., 2014b. “Numerical Analysis of Hydrocyclones with Different Conical Section Designs.” Minerals Engineering 62: 74–84. doi:10.1016/j.mineng.2013.12.003. Gibson, M. M., and Launder, B. E., 1978. “Ground Effects on Pressure Fluctuations in the Atmospheric Boundary Layer.” Journal of Fluid Mechanics 86 (3). Cambridge University Press: 491. doi:10.1017/S0022112078001251. Griffiths, W.D., and Boysan, F., 1996. “Computational Fluid Dynamics (CFD) and Empirical Modelling of the Performance of a Number of Cyclone Samplers.” Journal of Aerosol Science 27 (2). Pergamon: 281–304. doi:10.1016/0021-8502(95)00549-8. Habibian, M., Pazouki, M., Ghanaie, H., and Abbaspour-Sani, K., 2008. “Application of Hydrocyclone for Removal of Yeasts from Alcohol Fermentations Broth.” Chemical 155 Engineering Journal 138 (1–3): 30–34. doi:10.1016/j.cej.2007.05.025. He, P, Salcudean, M., and Gartshore, I. S., 1999. “A Numerical Simulation of Hydrocyclones.” Chemical Engineering Research and Design 77 (July): 429–41. doi:10.1205/026387699526412. Hirai, S., Takagi, T., and Matsumoto, M., 1988. “Predictions of the Laminarization Phenomena in an Axially Rotating Pipe Flow.” Journal of Fluids Engineering 110 (4). American Society of Mechanical Engineers: 424. doi:10.1115/1.3243573. Howard, J. H. G., Patankar, S. V., and Bordynuik, R. M., 1980. “Flow Prediction in Rotating Ducts Using Coriolis-Modified Turbulence Models.” Journal of Fluids Engineering 102 (4). American Society of Mechanical Engineers: 456. doi:10.1115/1.3240725. Hsieh, K. T., and Rajamani, R. T., 1991. “Mathematical Model of the Hydrocyclone Based on Physics of Fluid Flow.” AIChE Journal 37 (5): 735–46. Hsieh, KT. 1988. “Phenomenological Model of the Hydrocyclone.” University of Utah. Hwang, K.J., Hwang, Y.W. and Yoshida, H., 2013. “Design of Novel Hydrocyclone for Improving Fine Particle Separation Using Computational Fluid Dynamics.” Chemical Engineering Science 85: 62–68. doi:10.1016/j.ces.2011.12.046. Hwang, K.J., Wu, W.H., Qian, S., and Nagase, Y., 2008. “Separation Science and Technology CFD Study on the Effect of Hydrocyclone Structure on the Separation Efficiency of Fine Particles CFD Study on the Effect of Hydrocyclone Structure on the Separation Efficiency of Fine Particles.” Nagase Separation Science and Technology 4315: 3777–97. doi:10.1080/01496390802286637. Jones, W.P, and Launder, B. E., 1972. “The Prediction of Laminarization with a Two-Equation Model of Turbulence.” International Journal of Heat and Mass Transfer 15 (2): 301–14. doi:10.1016/0017-9310(72)90076-2. Kelsall, D.F. 1953. “A Further Study of the Hydraulic Cyclone.” Chemical Engineering Science 156 2 (6): 254–72. doi:10.1016/0009-2509(53)80044-8. Kelsall, DF. 1952. “A Study of the Motion of Solid Particles in a Hydraulic Cyclone.” Trans Inst Chem Eng 30: 87–108. Kikuyama, K., Murakami, M., and Nishibori, K., 1983. “Development of Three-Dimensional Turbulent Boundary Layer in an Axially Rotating Pipe.” Journal of Fluids Engineering 105 (2). American Society of Mechanical Engineers: 154. doi:10.1115/1.3240955. Kikuyama, K., Murakami, M., Nishibori, K., and Maeda, K., 1983. “Flow in an Axially Rotating Pipe : A Calculation of Flow in the Saturated Region.” Bulletin of JSME 26 (214). The Japan Society of Mechanical Engineers: 506–13. http://ci.nii.ac.jp/naid/110002364025/en. Knowles, S. R., Woods, D. R., Feuerestein, I. A., 1973. “The Velocity Distribution within a Hydrocyclone Operating Operating without an Air Core.” The Candian Journal of Chemical Engineering 51 (June): 263–71. Kolmogorov, A. N. 1941. “Local Structure of Turbulence in Incompressible Viscous Fluid Foir Very Large Reynolds Number.” Doklady Akademia Nauk SSSR 30 (299–303). Kuang, S. B., Chu, K. W., Yu, A. B., and Vince, A., 2012. “Numerical Study of Liquid-Gas-Solid Flow in Classifying Hydrocyclones: Effect of Feed Solids Concentration.” Minerals Engineering 31: 17–31. doi:10.1016/j.mineng.2012.01.003. Launder, B. E., Priddin, C. H., and Sharma, B. I., 1977. “The Calculation of Turbulent Boundary Layers on Spinning and Curved Surfaces.” Journal of Fluids Engineering 99 (1). American Society of Mechanical Engineers: 231. doi:10.1115/1.3448528. Launder, B. E., and Shima, N.. 1989. “Second-Moment Closure for the near-Wall Sublayer - Development and Application.” AIAA Journal 27 (10): 1319–25. doi:10.2514/3.10267. Launder, B.E. 1989. “Second-Moment Closure and Its Use in Modelling Turbulent Industrial Flows.” International Journal for Numerical Methods in Fluids 9 (8): 963–85. 157 Launder, B.E., and Sharma, B. I., 1974. “Application of the Energy-Dissipation Model of Turbulence to the Calculation of Flow near a Spinning Disc.” Letters in Heat and Mass Transfer 1 (2): 131–37. doi:10.1016/0094-4548(74)90150-7. Launder, B E, Reece, G. J. and Rodi, W., 1975. “Progress in the Development of a Reynolds-Stress Turbulence Closure.” J. Fluid Meeh 68 (3): 537–66. doi:10.1017/S0022112075001814. Launder, B. E. 1989. “Second-Moment Closure: Present… and Future?” International Journal of Heat and Fluid Flow 10 (4): 282–300. doi:10.1016/0142-727X(89)90017-9. Leonard, A. 1975. “Energy Cascade in Large-Eddy Simulations of Turbulent Fluid Flows.” Advances in Geophysics 18 (PART A): 237–48. doi:10.1016/S0065-2687(08)60464-1. Lien, F. S, and Leschziner, M. A., 1994. “Assessment of Turbulence-Transport Models Including Non-Linear Rng Eddy-Viscosity Formulation and Second-Moment Closure for Flow over a Backward-Facing Step.” Computers & Fluids 23 (8). Pergamon: 983–1004. doi:10.1016/0045-7930(94)90001-9. Lilly, D. K. 1992. “A Proposed Modification of the Germano Subgrid‐scale Closure Method.” Physics of Fluids A: Fluid Dynamics 4 (3). American Institute of PhysicsAIP: 633–35. doi:10.1063/1.858280. Ma, L., Ingham, D. B., and Wen, X., 2000. “Numerical Modelling of the Fluid and Particle Penetration through Small Sampling Cyclones.” Journal of Aerosol Science 31 (9): 1097–1119. doi:10.1016/S0021-8502(00)00016-1. Mackenzie, J. A., 2014. “The Effects of Turbulent Drag Reducing Additives on Hydrocyclone Operation.” University of British Columbia. MacKenzie, J., Martinez, D. M., and Olson, J. A.. 2014. “A Quantitative Analysis of Turbulent Drag Reduction in a Hydrocyclone.” Canadian Journal of Chemical Engineering 92 (8): 1432–43. doi:Doi 10.1002/Cjce.21989. 158 MacNaughton, J. 1906. No Title, issued 1906. Malhotra, A., Branion, R. M. R., and Hauptmann, E. G., 1994. “Modelling the Flow in a Hydrocyclone.” The Canadian Journal of Chemical Engineering 72 (6). Wiley Subscription Services, Inc., A Wiley Company: 953–60. doi:10.1002/cjce.5450720603. Mansour, N. N., Kim, J., and Moin, P., 1988. “Reynolds-Stress and Dissipation-Rate Budgets in a Turbulent Channel Flow.” Journal of Fluid Mechanics 194 (1). Cambridge University Press: 15. doi:10.1017/S0022112088002885. Mokni, I., Dhaouadi, H., Bournot, P., and Mhiri, H., 2015. “Numerical Investigation of the Effect of the Cylindrical Height on Separation Performances of Uniflow Hydrocyclone.” Chemical Engineering Science 122. Elsevier: 500–513. doi:10.1016/j.ces.2014.09.020. Monredon, T. C., Hsieh, K. T., and Rajamani, R. K., 1992. “Fluid Flow Model of the Hydrocyclone: An Investigation of Device Dimensions.” International Journal of Mineral Processing 35 (1–2): 65–83. doi:10.1016/0301-7516(92)90005-H. Mousavian, S. M., and Najafi, A. F., 2009. “Influence of Geometry on Separation Efficiency in a Hydrocyclone.” Arch Appl Mech 79: 1033–50. doi:10.1007/s00419-008-0268-8. Murakami, M., and Kikuyama, K., 1980. “Turbulent Flow in Axially Rotating Pipes.” Journal of Fluids Engineering 102 (1). American Society of Mechanical Engineers: 97. doi:10.1115/1.3240633. Murphy, S., Delfos, R., Pourquié, M. J. B. M., Olujić, Ž., Jansens, P. J., and Nieuwstadt, F.T.M., 2007. “Prediction of Strongly Swirling Flow within an Axial Hydrocyclone Using Two Commercial CFD Codes.” Chemical Engineering Science 62 (6): 1619–35. doi:10.1016/j.ces.2005.10.031. Murthy, Y. R., and Udaya Bhaskar, K., 2012. “Parametric CFD Studies on Hydrocyclone.” Powder Technology 230. Elsevier B.V.: 36–47. doi:10.1016/j.powtec.2012.06.048. Narasimha, M., Brennan, M., and Holtham, P. N., 2006. “Large Eddy Simulation of 159 Hydrocyclone—prediction of Air-Core Diameter and Shape.” International Journal of Mineral Processing 80 (1): 1–14. doi:10.1016/j.minpro.2006.01.003. Narasimha, M, Brennan, M., and Holtham, P. N., 2007. “A Review of CFD Modeling for Performance Predictions of Hydrocyclone.” Engineering Applications of Computational Fluid Mechanics 1 (2): 109–25. doi:10.1080/19942060.2007.11015186. Nishibori, K., Kikuyama, K., and Murakami, M., 1987. “Laminarization of Turbulent Flow in the Inlet Region of an Axially Rotating Pipe.” JSME International Journal 30 (260). The Japan Society of Mechanical Engineers: 255–62. doi:10.1299/jsme1987.30.255. Nowakowski, A. F., Cullivan, J. C., Williams, R. A., and Dyakowski, T., 2004. “Application of CFD to Modelling of the Flow in Hydrocyclones. Is This a Realizable Option or Still a Research Challenge?” Minerals Engineering 17 (5): 661–69. doi:10.1016/j.mineng.2004.01.018. Nowakowski, A. F., and Dyakowski, T., 2003. “Investigation of Swirling Flow Structure in Hydrocyclones.” Chemical Engineering Research and Design 81 (8): 862–73. doi:10.1205/026387603322482103. Ohasi, H., Maeda, S., 1958. “Motion of Water in a Hydraulic Cyclone.” Chemical Engineering Japan 22 (4): 200–207. Ohtake, T., Usuda, M., and Kadoya, T., 1987. “A Fundamental Study of Hydrocyclones Part 1. Flow Pattern in the Hydrocyclone.” Japan TAPPI Journal 41 (2): 60–64. Olson, J. A., and Kerekes, R. J., 2017. “The Motion of Fibres in Turbulent Flow.” J. Fluid Mech 377. Cambridge University Press: 47–64. doi:10.1017/S0022112098002973. Orlandi, P., and Fatica, M., 1997. “Direct Simulations of Turbulent Flow in a Pipe Rotating about Its Axis.” Journal of Fluid Mechanics 343 (July). Cambridge University Press: 43–72. doi:10.1017/S0022112097005715. Pericleous, K.A., and Rhodes, N., 1986. “The Hydrocyclone Classifier — A Numerical 160 Approach.” International Journal of Mineral Processing 17 (1). Elsevier: 23–43. doi:10.1016/0301-7516(86)90044-X. Petty, C. A., and Parks, S. M., 2004. “Flow Structures within Miniature Hydrocyclones.” Minerals Engineering 17 (5): 615–24. doi:10.1016/j.mineng.2004.01.020. Quian, L., Changlie, D., Jirun, X., Lixin, Y., and Guangai, X., 1989. “Comparison of the Performance of Water-Sealed and Commercial Hydrocyclones.” International Journal of Mineral Processing Elsevier Science Publishers B.V 25: 297–310. Reich, G, and Beer, H., 1989. “Fluid Flow and Heat Transfer in an Axially Rotating Pipe-L. Effect of Rotation on Turbulent Pipe Flow” 32 (3): 551–62. http://ac.els-cdn.com/0017931089901439/1-s2.0-0017931089901439-main.pdf?_tid=4e9ceb60-ffcc-11e6-94aa-00000aab0f02&acdnat=1488516543_287bc25f0458e8555db7a8d3cb610381. Renner, V. G., Cohen, H. E., 1978. “Measurement and Interpretation of Size Distribution of Particles within a Hydrocyclone.” TRANSACTIONS OF THE INSTITUTION OF MINING AND METALLURGY SECTION C-MINERAL PROCESSING AND EXTRACTIVE METALLURGY 87 (JUN): C139–45. Rietema, K. 1961a. “Performance and Design of hydrocyclones—III.” Chemical Engineering Science 15 (3–4): 310–19. doi:10.1016/0009-2509(61)85035-5. ———. 1961b. “Performance and Design of hydrocyclones—IV.” Chemical Engineering Science 15 (3–4): 298–302. doi:10.1016/0009-2509(61)85036-7. Rietema, K., 1961. “Performance and Design of hydrocyclones—II : Pressure Drop in the Hydrocyclone.” Chemical Engineering Science 15 (3). Pergamon: 303–9. doi:10.1016/0009-2509(61)85034-3. Rodi, W. 1980. “Turbulence Models and Their Application in Hydraulics – A State of the Art Review.” In IAHR. Delft, The Netherlands. Roldan-Villasana, E. J., Williams, R. A., and Dyakowski, T., 1993. “The Origin of the Fish-161 Hook Effect in Hydrocyclone Separators.” Powder Technology 77 (3): 243–50. doi:10.1016/0032-5910(93)85017-4. Rotta, J. 1951. “Statistical Theory of Nonhomogeneous Turbulence.” London. Saffman, P. G. 1974. “Model Equations for Turbulent Shear Flow.” Studies in Applied Mathematics 53 (1): 17–34. doi:10.1002/sapm197453117. Sagaut, Pierre. 2001. Large Eddy Simulation for Incompressible Flows. Scientific Computation. Berlin, Heidelberg: Springer Berlin Heidelberg. doi:10.1007/978-3-662-04416-2. Saidi, Maysam, Reza Maddahian, and Bijan Farhanieh. 2013. “Numerical Investigation of Cone Angle Effect on the Flow Field and Separation Efficiency of Deoiling Hydrocyclones.” Heat and Mass Transfer/Waerme- Und Stoffuebertragung 49 (2): 247–60. doi:10.1007/s00231-012-1085-8. Saidi, Maysam, Reza Maddahian, Bijan Farhanieh, and Hossein Afshin. 2012. “Modeling of Flow Field and Separation Efficiency of a Deoiling Hydrocyclone Using Large Eddy Simulation.” International Journal of Mineral Processing 112: 84–93. doi:10.1016/j.minpro.2012.06.002. Sevilla, E. M., and R. M. R. Branion. 1997. “The Fluid Dynamics of Hydrocyclones.” Journal of Pulp and Paper Science 23 ((2)): 85–93. Slack, M.D., R.O. Prasad, A. Bakker, and F. Boysan. 2000. “Advances in Cyclone Modelling Using Unstructured Grids.” Chemical Engineering Research and Design 78 (8). Elsevier: 1098–1104. doi:10.1205/026387600528373. Smagorinsky, J. 1963. “GENERAL CIRCULATION EXPERIMENTS WITH THE PRIMITIVE EQUATIONS.” Monthly Weather Review 91 (3): 99–164. doi:10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2. Speziale, Charles G., Sutanu Sarkar, and Thomas B. Gatski. 1991. “Modelling the Pressure–strain Correlation of Turbulence: An Invariant Dynamical Systems Approach.” Journal of 162 Fluid Mechanics 227 (1). Cambridge University Press: 245. doi:10.1017/S0022112091000101. Svarovsky, L. 1984. Hydrocyclones. Holt, Rinehart and Winston. Versteeg, H. K.;, and W.; Malalasekera. 2007. Introduction to Computational Fluid Dynamics. Pearson Education. 2nd ed. Vol. 2nd editio. Harlow: Pearson Education. doi:10.2514/1.22547. Vieira, L. G M, E. A. Barbosa, J. J R Damascene, and M. A S Barrozo. 2005. “Performance Analysis and Design of Filtering Hydrocyclones.” Brazilian Journal of Chemical Engineering 22 (1): 143–52. doi:10.1590/S0104-66322005000100015. Vieira, L. G M, Beatriz C. Silvério, J. J R Damasceno, and M. A S Barrozo. 2011. “Performance of Hydrocyclones with Different Geometries.” Canadian Journal of Chemical Engineering 89 (4): 655–62. doi:10.1002/cjce.20461. Wang, B., and A. B. Yu. 2006. “Numerical Study of Particle-Fluid Flow in Hydrocyclones with Different Body Dimensions.” Minerals Engineering 19 (10): 1022–33. doi:10.1016/j.mineng.2006.03.016. ———. 2008. “Numerical Study of the Gas-Liquid-Solid Flow in Hydrocyclones with Different Configuration of Vortex Finder.” Chemical Engineering Journal 135 (1–2): 33–42. doi:10.1016/j.cej.2007.04.009. White, A.; 1964. “Flow of a Fluid in an Axially Rotating Pipe.” Journal of Mechanical Engineering Science 6 (1): 47–52. Wilcox, David C. 2006. Turbulence Modeling for CFD. 3rd ed. La Canada, California: DCW Industries. Xu, Y., Xingfu S., Ze Sun, Bo Tang, Ping Li, and Jianguo Yu. 2013. “Numerical Investigation of the Effect of the Ratio of the Vortex-Finder Diameter to the Spigot Diameter on the Steady State of the Air Core in a Hydrocyclone.” Industrial and Engineering Chemistry Research 163 52 (15): 5470–78. doi:10.1021/ie302081v. Yoshioka, N.; Hotta, Y. 1955. “Liquid Cyclone as a Hydraulic Classifier.” Chemical Engineering Japan 19 (12): 633–41.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Hydrodynamics of the hydrocyclone flow field : effects...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Hydrodynamics of the hydrocyclone flow field : effects of turbulence modelling, geometrical and design… Zaman, Ehsan 2017
pdf
Page Metadata
Item Metadata
Title | Hydrodynamics of the hydrocyclone flow field : effects of turbulence modelling, geometrical and design parameters |
Creator |
Zaman, Ehsan |
Publisher | University of British Columbia |
Date Issued | 2017 |
Description | This thesis describes a numerical study of the effect of swirl on the flow hydrodynamics in hydrocyclones and rotating pipes. Kinetic (force balance) and kinematic (velocity field) analysis were performed to identify the relative importance of turbulence, convective, and pressure forces in these flows. We showed that convective and pressure forces in hydrocyclones are three orders of magnitude larger than the size of turbulent forces, and seven orders of magnitude greater than the viscous forces. We derived a new dimensionless parameter based on a ratio of Euler numbers that gave a unique relationship for pressure along the axis of a hydrocyclone for all the cone angles tested. This finding enables rigorous scaling of hydrocyclones of differing cone angles. A complementary numerical study of flow in rotating pipes was carried out to elucidate the relative importance of convection and turbulence. We identified a dimensionless rotation parameter that delineates the condition at which decreasing turbulence force equals increasing convective force as rotational speed increases. This dimensionless number establishes a criterion for knowing which forces are dominant, and thereby a rational basis for choosing CFD (Computational Fluid Dynamics) models that are both cost-effective and accurate. Lastly, we conducted a sensitivity study to determine the effects of varying values of constants in the quadratic pressure-strain (QPS) Reynolds Stress turbulence model (RSM) for CFD modelling of swirl in hydrocyclones. The results identified the changes necessary to attain accurate predictions, such as overcoming the under-prediction problems of RSM. In addition, these findings laid the groundwork for a swirl-specific turbulence model for hydrocyclones. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2017-10-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0357149 |
URI | http://hdl.handle.net/2429/63326 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2017-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2017_november_zaman_ehsan.pdf [ 8.16MB ]
- Metadata
- JSON: 24-1.0357149.json
- JSON-LD: 24-1.0357149-ld.json
- RDF/XML (Pretty): 24-1.0357149-rdf.xml
- RDF/JSON: 24-1.0357149-rdf.json
- Turtle: 24-1.0357149-turtle.txt
- N-Triples: 24-1.0357149-rdf-ntriples.txt
- Original Record: 24-1.0357149-source.json
- Full Text
- 24-1.0357149-fulltext.txt
- Citation
- 24-1.0357149.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0357149/manifest