NUMERICAL INVESTIGATION OF ROTATING STALL IN CENTRIFUGAL COMPRESSORS by Taher Mohammed Abdelhakim Fattouh Halawa B.Sc., Cairo University, 2005 M.Sc., Cairo University, 2008 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) September 2016 © Taher Mohammed Abdelhakim Fattouh Halawa, 2016 ii Abstract Rotating stall is a complicated phenomenon which may occur when a centrifugal compressor works at low flow rate condition close to the compressor surge limit. Rotating stall and surge may cause a large amount of damage to the compressor components due to instability and high pressure fluctuations which increase vibration of the compressor blades. This research introduces numerical simulations of the rotating stall phenomenon in a high speed centrifugal compressor with vaned and vaneless diffusers. The present study discusses the improvement of two important methods to control stall: the air injection method and the casing grooves treatment method. The method of air injection was optimized by monitoring the effects resulting from the variation of the angle of injection and the injection flow rate. The efficiency of the casing grooves method was related to the grooves' locations, cross section dimensions, and the number of grooves. The current research also presents estimations of the deformation of the compressor blades during the transition stage between stall and surge in order to determine whether the blades will deform safely within the tip clearance gap without touching the casing. Rotating stall simulations showed that the initiation of stall is due to the flow separation and the formation of the tip leakage flow. Simulations showed also that stall development appears at the vaneless area between the impeller and diffuser. Air injection simulations showed that stability is optimized when the air injection angle is approximately double the value of the diffuser inlet vane angle. Simulations showed that in order to achieve effective flow circulation, the aspect ratio of the casing groove cross section should be at least equal to one. Structural analysis showed that the most dangerous region at the deformed blade’s tip surface is located at the transition area between the axial and radial blade profile. Results also showed that the deformation of the blades at the tip region can be limited by changing the shaft bearing stiffness; this is useful to prevent damage to the compressor blades. iii Preface A version of section 5.1 in Chapter 5 was collected from the data presented in the following articles: Halawa, T., Alqaradawi, M., Badr, O., and Gadala, M. S., "Numerical Simulation and Control of Rotating Stall in a Transonic Centrifugal Compressor," Proc. ASME 2014 12th Biennial Conference on Engineering Systems Design and Analysis, American Society of Mechanical Engineers, pp. V002T011A024-V002T011A024. I generated the numerical model, analyzed results, and wrote the manuscript in consultation with M. S., Gadala. M., Alqaradawi (Assistant professor at the Mechanical and Industrial Engineering Department at Qatar University) and O., Badr (Professor and Head of the Mechanical Engineering Department at the British University in Egypt) enhanced the paper’s quality by suggesting the addition of specific numerical model validations. Halawa, T., Alqaradawi, M., Badr, O., and Gadala, M. S., "Numerical Investigation of Rotating Stall Characteristics and Active Stall Control in Centrifugal Compressors," Proc. ASME 2014 Power Conference, American Society of Mechanical Engineers, pp. V002T011A002-V002T011A002. I created the hypothesis, designed the numerical model, analyzed the model outputs, and wrote the manuscript in consultation with M.S., Gadala. M., Alqaradawi and O., Badr contributed their experience in the discussion of the numerical results. Halawa, T. and Gadala, M. S., “Rotating Stall simulation in Axial and Centrifugal Compressors," Accepted for publication in the International Journal for Computational Methods in Engineering Science and Mechanics, July 2016. I generated the numerical model, analyzed results, and wrote the manuscript in consultation with M. S., Gadala. iv A version of section 5.2 in Chapter 5 was taken from the following article: Halawa, T., Alqaradawi, M., Shahin, I., Gadala, M. S., and Badr, O., "Numerical investigation of rotating stall in centrifugal compressor with vaned and vaneless diffuser." Journal of Thermal Science 24.4 (2015): 323-333. I performed the simulations for the compressor model with vaned diffusers and participated inthe simulations of the compressor model with vaneless diffusers. I wrote the manuscript in consultation with M. S., Gadala. Also, I., Shahin (Postdoctoral Fellow at the Mechanical and Industrial Engineering Department at Qatar University) performed the numerical simulations for most of the analysis of the compressor model with vaneless diffusers and analyzed the data. M., Alqaradawi and O., Badr checked the concepts and results for accuracy. A version of section 5.3 in Chapter 5 was collected from the data published in the following articles: Halawa, T., Gadala, M. S., Alqaradawi, M., and Badr, O., 2015, "Optimization of the Efficiency of Stall Control Using Air Injection for Centrifugal Compressors," Journal of Engineering for Gas Turbines and Power, 137(7), p. 072604. Halawa, T., Gadala, M., Alqaradawi, M., and Badr, O., "Numerical Simulation of Stall Development into Surge and Stall Control Using Air Injection in Centrifugal Compressors," Proc. Proc. ASME Power 2014 Conference, Baltimore, Maryland, ASME paper No. POWER2014-32053. Halawa, T., Alqaradawi, M., Badr, O., and Gadala, M. S., "Numerical Investigation of Steady air Injection Flow to Control Rotating Stall in Centrifugal Compressors," Proc. ASME 2014 12th Biennial Conference on Engineering Systems Design and Analysis, American Society of Mechanical Engineers, pp. V002T011A025-V002T011A025. I created the numerical analysis and the optimization, and wrote the manuscript in consultation with M. S., Gadala. Also, M., Alqaradawi and O., Badr were involved in the discussion of the final results. v A version of section 5.4 in Chapter 5 was collected from the data presented in the following articles: Halawa, T., Alqaradawi, M., Badr, O., and Gadala, M. S., "Numerical Simulation of Using Combined Active and Passive Stall Control Techniques in Centrifugal Compressors," Proc. ASME 2014 International Mechanical Engineering Congress and Exposition, American Society of Mechanical Engineers, pp. V007T009A026-V007T009A026. Halawa, T., Gadala, M. S., Alqaradawi, M., and Badr, O., 2016, "Influence of Changing Casing Groove Parameters on the Performance of Centrifugal Compressors Near Stall Condition." Journal of Fluids Engineering 138.2 (2016): 021104. I generated the design, perormed the simulations, extracted results, and wrote the manuscript in consultation with M. S., Gadala. Some modifications for the numerical model design were suggested by M., Alqaradawi and O., Badr. vi Table of Contents Abstract ................................................................................................................................................. ii Preface ..................................................................................................................................................iii Table of Contents ................................................................................................................................ vi List of Tables ....................................................................................................................................... ix List of Figures ....................................................................................................................................... x List of Symbols .................................................................................................................................. xiv Acknowledgements ........................................................................................................................... xvi Dedication ......................................................................................................................................... xvii 1) CHAPTER 1: INTRODUCTION ............................................................................................................. 1 1.1 Compressor Map and the Safe Operating Zone ...................................................................... 1 1.2 Overview of the Rotating Stall Concept ................................................................................. 2 1.3 Surge and its Relation to Rotating Stall .................................................................................. 3 1.4 Rotating Stall Control ............................................................................................................. 4 1.5 Scope of the Work .................................................................................................................. 6 2) CHAPTER 2: LITERATURE REVIEW ................................................................................................... 8 2.1 Rotating Stall Characteristics ................................................................................................. 8 2.2 Compressor Safe Margin Enhancement Methods ................................................................ 11 2.2.1 Casing Treatment .......................................................................................................... 11 2.2.2 Air Injection .................................................................................................................. 14 2.3 Structural Analysis of Compressor Blades ........................................................................... 16 2.4 Synopsis of the Literature ................................................................................................... 17 3) CHAPTER 3: NUMERICAL METHODOLOGY ..................................................................................... 19 3.1 Turbulence Modeling ............................................................................................................ 20 3.2 Near-Wall Treatments for Turbulent Flow ........................................................................... 20 3.2.1 Wall Functions .................................................................................................................... 21 3.3 Dynamic Mesh Model .......................................................................................................... 23 3.3.1 Dynamic Mesh Concept ..................................................................................................... 23 3.3.2 Dynamic Mesh Conservation Equations ............................................................................ 25 3.4 Fluid-Structure Interface ...................................................................................................... 26 vii 4) CHAPTER 4: SELECTED CASE STUDY MODELING AND VALIDATION ............................................. 28 4.1 Description of the Case Study .............................................................................................. 28 4.2 Numerical Modeling of the Case Study ................................................................................ 29 4.2.1 CFD Models .................................................................................................................. 29 4.2.2 Structural Model ........................................................................................................... 31 4.2.3 Meshing Analysis .......................................................................................................... 32 4.3 Validation of the Numerical Model ...................................................................................... 36 5) CHAPTER 5: RESULTS AND DISCUSSION .......................................................................................... 44 5.1 Rotating Stall Simulation for the Compressor with a Vaned Diffuser ................................. 44 5.1.1 Overall View of Stall Characteristics ............................................................................ 44 5.1.2 Tip Leakage Flow Variation With Time ....................................................................... 46 5.1.3 Rotating Stall Effects on the Compressor Blades ......................................................... 49 5.1.4 The Most Affected Region by Rotating Stall ................................................................ 51 5.2 Rotating Stall Simulation for a Compressor with a Vaneless Diffuser and its Comparison to a Compressor with a Vaned Diffuser ............................................................................................... 54 5.2.1 Monitoring the Fluctuating Parameters During Surge .................................................. 56 5.2.2 Velocity Contours Comparison ..................................................................................... 60 5.2.3 Pressure Contours Comparison ..................................................................................... 62 5.2.4 Velocity Vectors Comparison ....................................................................................... 63 5.3 Air Injection Simulations ...................................................................................................... 64 5.3.1 Description of the Location and Orientation of the Injected Air .................................. 65 5.3.2 Clarification of the Air Injection Effect ........................................................................ 66 5.3.3 Comparison Between all Injection Cases ...................................................................... 67 5.3.4. Comparison Between the Best Injection Cases .................................................................. 72 5.3.5. Effect of Changing the Injection Angle in the Radial-Axial Plane ................................... 78 5.4 Casing Treatment Simulations .............................................................................................. 80 5.4.1 Simulation of Different Casing Treatment Configurations ........................................... 80 5.4.2 Simulation of the Casing Grooves with the Variation of their Parameters ................... 85 5.5 Structural Analysis of the NASA CC3 Compressor ............................................................. 98 5.5.1 Stress Variation Near Surge for the Controlled and the Uncontrolled Compressors .......... 98 5.5.2 Total Force and Moment Variations Near Surge for the Controlled and the Uncontrolled Compressors ............................................................................................................................... 100 viii 5.5.3 Deflection of the Blades During Stall for the Controlled and the Uncontrolled Compressors .................................................................................................................................................... 102 5.5.4 Effect of Variation in the Bearing Stiffness on the Deformation of the Blades Tip Surfaces During Stall ................................................................................................................................ 104 6) CHAPTER 6: CONCLUSIONS & FUTURE WORK .............................................................................. 108 6.1 Conclusions ....................................................................................................................... 108 6.1.1 Conclusions on the Simulations of Rotating Stall for the Compressor With Vaned Diffuser.. ..................................................................................................................................... 108 6.1.2 Notes on the Difference Between Compressors With Vaned and Vaneless Diffusers in the Stall Development ................................................................................................................ 109 6.1.3 Comments on the Air Injection Simulations ............................................................... 109 6.1.4 Comments on the Casing Treatment Simulations ....................................................... 111 6.1.5 Conclusions on the Structural Analysis Work ............................................................ 112 6.2 Future Work ....................................................................................................................... 113 REFERENCES ...................................................................................................................................... 114 APPENDIX ........................................................................................................................................... 119 A.1 The Conservation Equation of Transport and its Discretization ......................................... 119 A.2 Boussinesq Hypothesis ....................................................................................................... 119 A.3 Realizable ݇ െ ߝ Model Equations ..................................................................................... 120 A.4 Definitions of the Near-Wall Treatment Parameters .......................................................... 120 A.5 Equations for the Using of Wall Functions for Momentum ............................................... 121 A.6 Equations for the Using of Wall Functions for Energy ...................................................... 122 ix List of Tables Table 2-1 Literature summary ............................................................................................................ 18 Table 4-1 NASA CC3 compressor specifications ............................................................................... 29 Table 5-1 Stiffness values .................................................................................................................. 105 x List of Figures Figure 1-1 Schematic of the compressor pressure map [1]. ................................................................... 2 Figure 1-2 Flow rate variation for mild and deep surge cases [1]. ........................................................ 4 Figure 1-3 Casing treatment of a centrifugal compressor [21]. ............................................................. 5 Figure 1-4 Schematic of an air injection for a centrifugal compressor [22]. ......................................... 6 Figure 1-5 Clarification of the stall and surge region on the compressor pressure map ........................ 6 Figure 2-1 Schematic view of the stall inception conditions. ................................................................ 8 Figure 2-2 Influence of changing the clearance ratio on the onset of instability [4]. .......................... 10 Figure 2-3 Effect of changing the grooves’ number and location on the surge margin and efficiency losses [39]. ........................................................................................................................................... 12 Figure 2-4 Air injection system used for the NASA CC3 compressor [22]. ....................................... 14 Figure 3-1 Overview of the density based method. ............................................................................. 19 Figure 3-2 Classification of the near-wall region. ............................................................................... 21 Figure 3-3 Schematic drawing for the various zones result due to a non-periodic interaction. ........... 24 Figure 3-4 Faces and zones created by the interaction between two cell zones. ................................. 25 Figure 4-1 NASA CC3 centrifugal compressor [74] ........................................................................... 28 Figure 4-2 The single passage numerical model .................................................................................. 30 Figure 4-3 The full annulus numerical model ...................................................................................... 30 Figure 4-4 The structural model .......................................................................................................... 31 Figure 4-5 Global mapping method [75] ............................................................................................. 32 Figure 4-6 Pressure variation at the mid span with the axial direction (z) .......................................... 33 Figure 4-7 Mass flow rate convergence for different CFL values ....................................................... 34 Figure 4-8 Meshing at the vaneless region ......................................................................................... 34 Figure 4-9 Meshing of the CFD model ............................................................................................... 35 Figure 4-10 Meshing of the mechanical model .................................................................................. 35 Figure 4-11 Compressor pressure and efficiency maps (numerical results versus measurements [69]).............................................................................................................................................................. 36 Figure 4-12 Pressure and Mach number contours for the current model results and for previous CFD model results [75] ................................................................................................................................ 37 Figure 4-13 Mean velocity variation with the % meridional distance (numerical results versus measurements [73]) .............................................................................................................................. 38 Figure 4-14 Velocity profile from hub to tip at the vaneless region at radius ratio of 1.1 (numerical results versus measurements [69]) ....................................................................................................... 38 Figure 4-15 Normalized pressure variation with the normalized meridional distance at the impeller tip (numerical results versus measurements [77]) ................................................................................ 39 Figure 4-16 Velocity at the diffuser different times of the surge cycle (numerical results versus measurements [72] ............................................................................................................................... 40 Figure 4-17 Stress distribution at impeller blade surface (The presented numerical model results compared with the NASA numerical model results [73]) .................................................................... 42 Figure 4-18 Axial deflection (for the presented numerical model results) .......................................... 43 xi Figure 5-1 Velocity vectors near impeller shroud during stall ............................................................ 44 Figure 5-2 Velocity vectors showing tip leakage flow ........................................................................ 45 Figure 5-3 Velocity vectors showing back flow impingement at impeller exit during stall ................ 45 Figure 5-4 Velocity contours near diffuser shroud surface during stall ............................................... 46 Figure 5-5 Path lines near shroud surface at different time frames ..................................................... 47 Figure 5-6 Velocity vectors at different planes through the impeller at an instant of time during stall.............................................................................................................................................................. 48 Figure 5-7 Pressure development during stall ...................................................................................... 49 Figure 5-8 Stream lines distribution during stall and surge at location near shroud surface ............... 50 Figure 5-9 Variation of the force affecting impeller blades with time during surge ........................... 51 Figure 5-10 Location of the numerical monitor points at the tip clearance gap .................................. 52 Figure 5-11 Dynamic pressure variation with time at points in the tip clearance gap ......................... 52 Figure 5-12 Pressure variation with time during stall and surge for a point in the vaneless region .... 53 Figure 5-13 Pressure monitoring with time at the vaneless region ...................................................... 54 Figure 5-14 The numerical model for the compressor with a vaneless diffuser .................................. 55 Figure 5-15 Description of the operating conditions used for the simulations of the vaneless and the vaned diffuser cases ............................................................................................................................. 56 Figure 5-16 Area weighted average of the mass flow rate at compressor exit near surge condition for the vaneless and vaned diffuser cases .................................................................................................. 57 Figure 5-17 Static pressure at radius ratio of 1.1 near surge condition for the vaneless and vaned diffuser cases ........................................................................................................................................ 58 Figure 5-18 Pressure surge cycle for the vaned diffuser case according to measurements [72] ......... 59 Figure 5-19 Dynamic pressure at radius ratio of 1.1 near surge condition for the vaneless and vaned diffuser cases ........................................................................................................................................ 59 Figure 5-20 Spectral analysis of the convergence history of the dynamic pressure at radius ratio of 1.1 near surge condition for the vaneless and vaned diffuser cases ..................................................... 60 Figure 5-21 Axial velocity contours upstream the impeller inlet at different operating conditions for the vaneless and the vaned diffuser cases ............................................................................................ 61 Figure 5-22 RMS of the static pressure at the impeller hub and blades surfaces at different operating conditions for the vaneless and the vaned diffuser cases. .................................................................... 62 Figure 5-23 Velocity vectors at 98% of span for the vaneless and vaned diffuser cases during surge 64 Figure 5-24 Location of the point selected for the air injection simulations ....................................... 65 Figure 5-25 Description of the injection locations and the injection angle ......................................... 66 Figure 5-26 Pressure contours near the shroud surface ....................................................................... 67 Figure 5-27 Velocity vectors near the shroud surface ......................................................................... 67 Figure 5-28 The number of diffuser passages with reversed flow for various angles of injection and injection mass flow rates ...................................................................................................................... 68 Figure 5-29 Flow path lines colored by Mach number value for various angles of injection at the minimum and maximum injection mass flow rate used ...................................................................... 69 Figure 5-30 Representation of one of the paths through a diffuser passage used to show results ....... 70 Figure 5-31 Mach number variation with the radial distance at different injection mass flow rates for injection angles of 10◦, 20◦, 30◦, 40◦ and 180◦ ....................................................................................... 71 xii Figure 5-32 Mach number and velocity angle variations for the best injection cases ........................ 73 Figure 5-33 Pressure contours (in Pascals) at impeller shroud for the best injection cases ................ 75 Figure 5-34 The CFD numerical model of the NASA CC3 compressor ............................................. 76 Figure 5-35 Variation of the force acting on blades with time for the best injection cases and for the case without using injection ................................................................................................................. 77 Figure 5-36 The angle of injection α at the radial-axial plane ............................................................. 78 Figure 5-37 Pressure variation at points in the vaneless region with time at Ɵ = 30◦ for α = 10◦, 15◦, 20◦ ......................................................................................................................................................... 79 Figure 5-38 Numerical models for three different casing treatment cases and for the case without casing treatment ................................................................................................................................... 81 Figure 5-39 Velocity contours at 0.99 of span for the four different cases ......................................... 82 Figure 5-40 Average Mach number variation with time at the position of the first slot ...................... 84 Figure 5-41 Average Mach number variation with time ...................................................................... 84 Figure 5-42 Isentropic total to total efficiency variation with time during stall for all cases .............. 85 Figure 5-43 The numerical model used for the simulations of the groove aspect ratio variation ........ 86 Figure 5-44 The numerical model used for the simulations of the groove location variation ............. 87 Figure 5-45 The numerical models used for the simulations of changing the number of grooves ...... 87 Figure 5-46 Velocity contours at a plane perpendicular to the axial direction with the clarification of the velocity vectors inside the groove for various groove aspect ratio values ..................................... 89 Figure 5-47 Radial velocity variation with the circumferential angle for 3 different groove aspect ratio values ........................................................................................................................................... 90 Figure 5-48 Mach number values at 99% of the impeller span for 3 different groove locations ........ 91 Figure 5-49 Static pressure variation in the circumferential direction at the blades leading edge at 99% of span .......................................................................................................................................... 92 Figure 5-50 Velocity stream lines at 98% of span for cases with different numbers of grooves ........ 93 Figure 5-51 Average velocity contours at the impeller meridional plane for cases with different numbers of grooves .............................................................................................................................. 94 Figure 5-52 Pressure variation with the normalized meridional distance for the cases of smooth casing and for the cases with 5 and 7 grooves ..................................................................................... 95 Figure 5-53 Total to total isentropic efficiency for cases with different numbers of grooves ............. 96 Figure 5-54 Velocity vectors inside the grooves and at a plane near the shroud surface at the best efficiency operating condition ............................................................................................................. 96 Figure 5-55 Velocity vectors inside the grooves and at a plane near the shroud surface at the stall operating condition .............................................................................................................................. 97 Figure 5-56 The variation of the peak value of the maximum principal stress with time for the impeller blade ....................................................................................................................................... 99 Figure 5-57 The variation of the peak value of equivalent stress with time for the impeller blade .. 100 Figure 5-58 The variation of forces and moments with time for both of the controlled and the uncontrolled cases .............................................................................................................................. 101 Figure 5-59 The maximum full tip blade surface deformation profile for the case of an uncontrolled compressor ......................................................................................................................................... 103 xiii Figure 5-60 The maximum full tip blade surface deformation profile for the case of controlled compressor using air injection at zone 2 ............................................................................................ 104 Figure 5-61 The location of the compressor shaft bearings and their numerical representation ....... 105 Figure 5-62 Deformation of the blade tip surface at zone 2 for different stiffness values ................ 106 Figure 5-63 Deformation of the blade tip surface at a specific part of zone 2 for different stiffness values ................................................................................................................................................. 107 xiv List of Symbols A Area vector [m2] fA Area vector at the face f [m2] jA Area vector for the face j pc Specific heat of fluid [J/kgK] C Turbulence model parameter E Empirical constant F Force [N] kG Kinetic energy production rate H Groove Height [m] K Turbulence kinetic energy per unit mass [m2/ s2] pk Turbulence kinetic energy at point P[m2/ s2] Ks Spring stiffness [N/m] M Moment [N.m] ṁ Mass flow rate [kg/s] nf Number of faces on the control volume Pr Molecular Prandtl number Prt Turbulent Prandtl number q Wall heat flux [Watt/ m2] R Radial distance [m] S Shear stress [N/m2] S Average shear stress [N/ m2] S The source term of T Time [s] PT Temperature at the cell adjacent to wall [K] wT Temperature at the wall [K] t Time step [s] t Incremental change in time [s] u Flow velocity vector [m/s] gu Mesh velocity of the moving mesh [m/s] ܷା Dimensionless velocity ܷ∗ Dimensionless velocity for CFD cU Mean velocity magnitude when * *Ty y pU Mean velocity of the fluid at point P [m/s] ఛܷ Frictional velocity[m/s] v Velocity vector xv V Mesh volume [m3] V Incremental change in mesh volume [m3] W Groove Width [m] X x- direction Y y-direaction py Distance from point P to the wall [m] ݕା Dimensionless distance from the wall *y Dimensionless distance from the wall for CFD *Ty Thermal conduction layer thickness Z Axial direction R Radial distance [m] Abbreviations CFD Computational fluid dynamics CFL Courant–Friedrichs–Lewy number FEM Finite Element Model FSI Fluid-Structure Interaction GCI Grid Convergence Index KSI Kilo pound per Square Inch NASA National Aeronautics and Space Administration NSE Navier Stokes Equations RMS Root Mean Square rpm Revolutions per minute SRCT Self-Recirculation Casing Treatment Greek letters ∅ Symbol for general flow variable Ε Turbulence dissipation rate [m2/ s3] Fluid density [Kg/m3] Diffusion coefficient The gradient of a value Ɵ Injection angle in the radial-tangential plane Α Injection angle in the radial-axial plane Von Kármán constant Dynamic viscosity [pa.s] t Turbulent viscosity [pa.s] the mean rate of rotation tensor xvi Acknowledgements In the name of Allah, I started this research hoping that it accelerates the wheel of progress in this field. I would like to express my sincere appreciation and infinite thanks to my supervisor, Professor Mohamed S. Gadala, who gave me a great scientific support and guidance to make progress in my research. His vast experience and unlimited co-operation and efforts have helped me to produce this research in the best way. Special thanks to the Qatar National Research Fund (QNRF) for supporting this research from the first moment of work to the end. I thank all members of my family, especially my wife, who helped me greatly through all of it, and my parents for their magnificent role of supporting and encouraging me to move ahead towards progress and success. xvii Dedication To my wife, my parents, and my children 1 1)CHAPTER 1: INTRODUCTION Centrifugal compressors are widely used in many applications, as turbochargers in automotive engines and gas turbine engines, in the processing of natural gas, and in refrigeration systems. The stable operating range of centrifugal compressors is wider than that of the axial compressors and the efficiency drop at low flow rates is lower. Also, the pressure ratio for a single stage is higher in centrifugal compressors compared to axial compressors, making centrifugal compressors more common in industrial applications these days. Over the last decade, several studies have discussed various ways to increase the stable working zone of centrifugal compressors by studying the reasons for the evolution of the aerodynamic instabilities represented by stall and surge, and developing modern control methods to detect these instabilities. 1.1 Compressor Map and the Safe Operating Zone Maximum efficiency can be achieved if the compressor operates at the design mass flow rate. If the compressor works at the off design condition, there is a possibility of the occurrence of instabilities. For any compressor, there is a specific stable region that extends from the surge limit at relatively low flow rates to the choke limit at relatively high flow rates. Near the surge limit zone, aerodynamic instabilities and rotating stall may occur and cause repetitive pressure pulses on the rotor blades, and this may result in damage to the blades. On the other hand, near the choke limit zone, there is a high probability of the sudden appearance of strong shock waves which can affect the blades. Figure 1-1 shows a schematic of the compressor pressure map with an indication of the compressor's stable margin, which is limited by the surge and choke lines. Most of the current research about rotating stall aims at increasing the compressor’s stable zone of operation by shifting the surge limit to the left as shown in Fig. 1-1. Some pressure maps refer to a line called the surge avoidance line, and this line is parallel to the surge limit inside the stable zone. The purpose of the surge avoidance line is to indicate the dangerous unstable zone close to the surge limit. 2 Figure 1-1 Schematic of the compressor pressure map [1]. 1.2 Overview of the Rotating Stall Concept The rotating stall is a flow disturbance where waves propagate at the rotor tip area at a speed of 50% to 75% of the rotor speed [2]. Many studies suggest that the early separation and the tip leakage flows are important reasons for stall inception [3-5]. The tip leakage flow is considered a secondary flow which separates from the main flow and moves through the tip clearance gap randomly, causing partial blocking of the impeller passages. The stall cells form due to the occurrence of the low velocity regions caused by the flow separation and the tip leakage flow effect. The rotating stall cells move in the same direction of the rotor blade motion corresponding to the inertial coordinates. The mechanism of stall inception is complicated and there are two main types of stall inception known as modes and spikes. Modes are flow disturbances formed near the blade tips with long length scales and small amplitudes, and move in the axial flow direction. Spikes are finite amplitude disturbances that propagate mainly in the rotor circumference. There is a considerable difference between modes and spikes in terms of their appearance time before stall initiation. Spikes can cause a stall within a few rotor revolutions but modes need hundreds of rotor revolutions to trigger a stall [6-8]. The slope of the pressure ratio curve with the flow rate 3 can indicate which type of stall inception is taking place; the curve slope is zero in the case of modes, while its value is negative in the case of spikes [9, 10]. Rotating stall may cause large losses in efficiency and compression performance. The existence of any stall cell may result in a sudden drop in the compressor’s total pressure ratio. Also, the formation of multiple stall cells minimizes the pressure ratio gradually. During stall, the average mass flow rate does not change but there are pressure fluctuations that affect the blades. As the stall develops with time, the amplitude of the pressure fluctuations increases and the temperature exceeds the designed value; this may cause fatal damage to the rotor blades. 1.3 Surge and its Relation to Rotating Stall As the flow separation and stalled area expand with time, surge may take place. Surge can be defined as strong fluctuations in the mass flow rate accompanied with a rise in temperature and pressure that may cause damage to the rotor seals and bearings and also to the rotor blades [11]. The main difference between stall and surge is that the average flow rate is not constant and reversed flow can occur during surge. The surge can increase the blades’ strain to the limit that the impeller blades may touch the casing surface or the stator vanes, causing significant damage to the compression unit [12]. Surge is classified into four different types depending on the level of pressure fluctuations [11]. The surge classes are mild, classic, modified and deep. For the mild surge, there are small pressure fluctuations without any back flow. The classic surge is similar to the mild surge in that there is no reversed flow but the amplitude of the pressure fluctuations is larger and occurs at lower frequencies compared to the mild surge. A more developed situation where the rotating stall combines with the classic surge, and the flow at the blades’ tip starts to move in the axial direction, is called a modified surge stage. If the pressure oscillations become larger and reversed flow forms, this can be described as deep surge. Figure 1-2 shows the difference between mild surge and deep surge in terms of the mass flow rate fluctuation overtime. It may be seen in Fig. 1-2 that the average mass flow rate is constant for the mild surge, but for the deep surge, the average mass flow rate is variable and the amplitude of fluctuations increases with time until the occurrence of the reversed flow. The rotating stall is often a precursor of a surge and the rotating stall may lead to a mild surge, which may then develop into a classic or deep surge. 4 (a) Mild surge (b) Deep surge Figure 1-2 Flow rate variation for mild and deep surge cases [1]. 1.4 Rotating Stall Control There are two main ways to control stall: the closed loop method and the open loop method. In the closed loop control method, there are some devices placed and fixed on the compressor casing in order to detect the rotating stall by measuring some of the flow parameters. The signals that come out of these devices are passed for filtering and then sent to a control unit in order to prevent or delay the onset of stall. In the open loop control method, the design of the compressor geometry is altered for the purpose of expanding the compressor’s stable zone by shifting the surge line. The open loop control method involves adjusting the inlet guide vanes, casing treatment and air injection control. The inlet guide vanes are angle-adjusted vanes which are placed before the impeller blades in order to control the flow angle at the impeller inlet when the flow rate is reduced [13, 14]. The casing treatment concept was first proposed by Hartmann et al. [15]. The idea of casing treatment control is to make a specific modification in the casing aimed at minimizing the tip leakage flow and increasing the kinetic energy in the stalled regions depending on the circulation effect. The most common types of casing slots are the casing discrete slots which are distributed in the circumferential direction and the continuous circumferential grooves which are designed to have the same cross section circumferentially [16, 17]. Figure 1-3 describes the casing discrete slots type for a turbocharger compressor. It may be noted from Fig. 1-3 that there 5 are some slots made in the casing near the impeller inlet region and the function of these slots is to transfer the tip leakage flow from the impeller tip area to the impeller inlet at a low flow rate operating condition. One disadvantage with the casing treatment is that the flow movement inside slots or grooves causes a drop in the efficiency due to the increased frictional interaction [18-20]. Figure 1-3 Casing treatment of a centrifugal compressor [21]. The air injection is used to decrease the flow incidence angle, and adding kinetic energy when the operating condition is corresponding to low mass flow rate. In the air injection method, pressurized air is injected through injectors mounted on the compressor casing at the impeller inlet or at the vaneless region. Figure 1-4 shows the air injection system of a centrifugal compressor from a previous experimental study [22]. The air is injected through eight nozzles fixed at the shroud surface of the impeller at the vaneless region area in an opposite direction to the impeller rotation direction (Fig. 1-4). The injection angle and mass flow rate can be adjusted by these injectors in order to enhance the injection efficiency. The air injection method has proved successful in decreasing the surge zone for both axial and centrifugal compressors [22-24]. Casing slots 6 Figure 1-4 Schematic of an air injection for a centrifugal compressor [22]. 1.5 Scope of the Work The scope of the work in this thesis is highlighted in the following objectives: Accurate simulation of the rotating stall phenomenon for a high speed centrifugal compressor and providing detailed explanations about the stall initiation and its development with time. Increasing the compressor stable margin by making an optimization to the factors which guide the air injection control method and the casing treatment control method (Fig. 1-5). Decreasing this region is required in order to increase the compressor stable zone. Figure 1-5 Clarification of the stall and surge region on the compressor pressure map 7 Expanding the information provided by previous studies regarding the structural analysis of the selected compressor by: Studying the dynamic structural behavior of the compressor during stall/surge by considering one way Fluid Structure Interaction (FSI) analysis. Assessment of blade deflection and stresses during stall/surge & its implication on the design stage. Finally, this research aims to provide a comparison between the behavior of controlled and uncontrolled compressors regarding the maximum deformation at the blades tips when using various values of the compressor shaft bearing stiffness based on the results of the dynamic structural model. In order to attain the above objectives, the following tasks are considered: i) Construction of CFD model for a high speed centrifugal compressor that is capable of simulating the compressor’s unstable behavior at high pressures and at low mass flow rates. ii) Developing the computational mesh for a section model of the compressor containing one blade passage, and testing the mesh density’s effect on the solution convergence in order to optimize the model accuracy. iii) Building the full annulus model containing all impeller blades and diffuser vanes, obtaining results concerning the identification of rotating stall parameters, and recording the early signs of stall formation. iv) Simulating modern active control methods and optimizing control parameters in order to find the best suitable combination of control actions to delay or prevent stall. v) Performing structural analysis on the deflection of impeller blades during stall by analyzing the FEM model results, which are obtained after transferring the pressure loads from the CFD model. 8 2)CHAPTER 2: LITERATURE REVIEW This chapter presents a survey of the literature for three main subjects: the characteristics of rotating stall, the methods used to increase the compressor's stable margin, and the structural analysis of the compressor blades. 2.1 Rotating Stall Characteristics The complexity of the rotating stall phenomenon led to much research about the reasons for stall initiation and the stall inception mechanism. Vo et al. [25] explained the mechanism of stall inception and the main reasons for stall initiation in compressors. It was found that flow separation causes the formation of secondary flow which moves between any two adjacent impeller passages through the casing clearance; this kind of flow is called the tip clearance flow. This study also clarified that part of the flow at the impeller exit turns back and hits the impeller blades randomly at the beginning of stall. It was found that the rotating stall starts when the tip leakage flow appears at the same time as the impeller blades are impinged with some back flow. Figure 2-1 indicates the two conditions which cause stall. Figure 2-1 Schematic view of the stall inception conditions. 9 The stall inception mechanism clarified by Vo et al. [25] is one of the most famous explanations of stall initiation in axial compressors. The stall inception mechanism in centrifugal compressors is close to but not identical to that for axial compressors; it has a more complex nature due to the flow direction at the compressor exit compared to the inlet. There is no complete explanation for the stall development criteria in centrifugal compressors. One of the main aims of the research presented in this thesis is to simulate and analyze the flow variations during stall, as a step towards improving some deficiencies in the understanding of stall complexity in centrifugal compressors. Pullan et al. [26] numerically and experimentally investigated stall inception caused by spikes, and examined their structure. It was proven that the starting point of stall cell formation is at the impeller leading edge, due to the effect of the flow with a high incedence angle on the vortices' generation rate. Results showed that there is a specified path for most of the vortices which moves from the suction side to the pressure side of the blade gradually, towards the shroud surface. Several studies related the stall progress rate to the flow at the blades tips where the stall cells form and develop with time. Huang et al. [27] created numerical simulations of the unsteady characteristics of rotating stall for a high speed centrifugal compressor. Numerical results showed that the impeller tip exit region is the location that is most affected by the pressure variation during stall. By monitoring pressure fluctuations at the blades tip area, it was found that frequency of the tip leakage vortex is close to half of the impeller blade passing frequency value. Geng et al. [28] studied the variation of the tip leakage flow overtime for an axial compressor. Results showed that the tip clearance flow propagates circumferencially with time during stall, and blocks some of the impeller passages. It was also found that the main flow passing through the impeller is affected by a reversed flow which is separated from the tip leakage flow during stall development. Yamada et al. [29] showed that the flow separation at the impeller leading edge generates votices which have a tornado shape, and confirmed that these vortices are the reason for the tip leakage flow formation at the beginning of a stall. Tomita et al. [30] presented a numerical and experimental investigation of the tip leakage vortex formation and its effect on stability. It was found that the link between stall and surge depends on the tip leakage vortex and its effect on limiting the safe operating range. 10 The rotating stall effects may be minimized by modifying or enhancing the design parameters of the compressor such as the clerance gap size, blade shape and blade angles. Schleer and Abhari [4] demonstrated the influence of modifying the clerance gap dimension between the blade tips and casing on the performance of a centrifugal compressor at stall condition. Results indicated that the compressor stable zone can be increased and the blade loading can be decreased by minimizing the tip clearance gap. The effect of the clearance ratio variation on the onset of instability is clarified on the pressure map shown in Fig. 2-2. The clearance ratio is defined as the ratio between the tip clearance value at impeller exit, and the diffuser height. Figure 2-2 Influence of changing the clearance ratio on the onset of instability [4]. Ramakrishna and Govardhan [31] compared the performance of an axial compressor with forward swept impeller blades with another compressor with unswept blades when the operating condition is close to stall. Based on the comparison between the two compressors, it was found that the diffusion factor at the impeller exit has lower values for the compressor with swept blades. It was also shown that the efficiency losses are lower for an increased clearance gap for swept impeller blades. Results confirmed that stall impact can be minimized by the sweeping of the impeller blades. Vagani et al. [32] performed various numerical simulations of the rotating stall phenomenon for centrifugal compressors and compared the numerical results with previously done meausrements. The numerical models covered different designs for the diffuser 11 vanes including the vanes’ shape and length. Lastly, it was shown that the numerical models provided accurate results comparedwith measurements for all cases. The work done by Schleer and Abhari [4], and Ramakrishna and Govardhan [31] prove that the tip leakage flow can be minimized by adjusting the tip clearance value and by sweeping the compressor blades, but there is a negative point where the efficiency losses are relatively high. For this reason, the current research presented in this thesis focuses on other methods of controlling the tip leakage flow and discusses ways to improve these methods as discussed in section 2.2; these methods provide relatively lower efficiency losses. 2.2 Compressor Safe Margin Enhancement Methods There are two significant methods used for increasing the safe margin of the compressor: the casing treatment method, and the air injection method. In the casing treatment method, the compressor casing design is modified by adding slots or grooves in order to circulate the weak flow near the blades tip. In the air injection method, the air is injected with a certain mass flow rate at specific locations from the shroud or the hub surfaces for the purpose of increasing the kinetic energy and modifying the flow angle at the stagnant areas. 2.2.1 Casing Treatment The first common type of casing treatment is called casing grooves, for which circumferentially continuous grooves are cut into the end wall at the shroud surface. The casing grooves lead to lower efficiency losses compared to other casing treatment types [33-35]. Kurokawa et al. [36] introduced an experimental and numerical study about using radial casing grooves in the compressor vaneless area. Results of this study indicated that the circulation caused by the casing grooves at the diffuser vaneless area modified the flow angle at the impeller exit during stall, and corrected the deviated flow direction. Back flow from the diffuser was also limited as a result of using grooves upstream of the diffuser entrance. Park et al. [19] presented a numerical model for the ring groove design in a centrifugal compressor. The results of this numerical model showed that the ring groove casing treatment increased the surge margin by 6% and caused the efficiency to drop by 6% relative to a compressor with smooth casing. Kim et al. [37] performed numerical simulations to test the optimal dimensions of the casing grooves and the tip clearance value in order to increase the stable opearting range. The 12 optimized model achieved an increase of about 5% in the surge margin but it was accompanied with 0.18% loss in efficiency relative to another model without casing treatment. Xi et al. [38] created a numerical optimization for the casing groove location at the diffuser of a centrifugal compressor. It was found that placing the casing grooves at the area that extends from the middle of the diffuser length to the diffuser exit can minimize stall in a more effective way than placing the grooves at the diffuser inlet. Ivor Day [39] constructed an experimental investigation to identify the best location for the casing grooves in order to delay the stall inception for an axial compressor. The mid of the chord distance was found to be the most successful position to place the grooves in terms of shifting the stall onset. Figure 2-3 indicates the impact of varying the location and number of grooves on the gain in surge margin and the losses in efficiency. This study found that the summation of the losses in efficiency for single groove components is larger than the total efficiency losses when using these components together as multiple grooves at the same time. Figure 2-3 Effect of changing the grooves’ number and location on the surge margin and efficiency losses [39]. To conclude, it is important to notice that these studies concerning the use of casing grooves did not consider the effect of combining all the factors of casing groove redesign. For 13 example, the work of Kim et al. [37] only considered the factor of changing the casing grooves’ dimensions, while the study of Xi et al. [38] mentioned the optimization of the casing grooves’ location. The impact of changing the number of casing grooves, their location and their dimensions is taken into consideration in this presented thesis. The second type of casing treatment is called discrete slots. For this type, a number of circumferentially or axially discrete slots are cut into the end wall at the shroud surface [40]. The efficiency losses caused by the discrete slots are greater than those for the casing grooves [41]. Lu et al. [42] numerically investigated the flow behaviorin a subsonic compressor enhanced by axial skewed slots. Results showed that the axial skewed slots increased the safe margin by shifting the vortices, caused by the tip leakage flow, downstream towards the rotor exit. Hembera et al. [18] simulated the flow inside a specific shape of axial slots for a transonic axial compressor. These slots have a semicircular shape and were set in the area between the inlet guide vanes and impeller at various angles to the axial direction. The circulated flow strength inside the semicircular slots was found to be stronger than for the regular slot shape. It was noted that there is a steady vortex motion inside the casing slots for the suction side of the inlet guide vanes close to their trailing edge. There is a special type of casing treatment that is called the holed casing treatment. The idea of the holed casing treatment is to make a number of cylindrical cavities in the casing connecting the impeller inlet to specific locations at the impeller shroud. The tip leakage flow passes through the holes in the casing and minimizes stall cells [20]. The holes’ diameter is an important parameter and it is related to the flow rate and flow direction inside the casing cavities [43]. Ding et al. [44] experimentally investigated the flow in a centrifugal compressor with the holed casing treatment. Measurements showed that there is a benefit to using the holed casing treatment design at both the low and high flow rates. At the high flow rate condition, there is a bypass flow from the impeller inlet to the impeller tip area while at the low flow rate condition, the tip leakage flow is drawn from the bleeding points at the impeller shroud to the impeller inlet. This flow mechanism was proven to increase the compressor’s stable region with relatively small efficiency losses. The Self Recirculation Casing Treatment (SRCT) method was used also recently and achieved success in increasing the compressor’s stability during stall, like the casing grooves and axial skewed slots [45]. 14 It is important to mention that the holed casing treatment method as described by Ding et al. [44] and the self recirculation casing treatment method as presented by Yang et al. [45], have practical manufacturing difficulties compared to other methods. For this reason, the casing grooves method was chosen to be studied in the current thesis instead of the other mentioned methods. This thesis offers a detailed discussion of the ways to improve the performance of the casing grooves method. 2.2.2 Air Injection Skoch [22] performed a series of experiments on the NASA CC3 centrifugal compressor equipped with an air injection system. The pressurized air was injected through nozzles which are fixed at the shroud surface between the impeller and diffuser. The tested range of the injection flow rate was 0.5% to 4% of the compressor inlet flow rate. The injection stream was directed circumferentially in the same and opposite directions of the impeller rotation. Figure 2-4 shows the air injection system of the NASA CC3 compressor. (a) Cross sectional view of the compressor indicating the injector (b) Schematic view for the injection orientations Figure 2-4 Air injection system used for the NASA CC3 compressor [22]. 15 It was found that the surge margin can be increased to 13.4% with an injection to inlet flow rate ratio of 0.9% if the injection has the same direction of impeller rotation. Experimental results showed that the efficiency losses are greater in the case of injecting the air in the reverse direction of impeller rotation. Further experiments done on the NASA CC3 compressor revealed that injecting air at the hub surface is not as effective as for the shroud surface injection, and that the injection towards the pressure side of the diffuser vanes is the best orientation that can maximize the surge margin [46]. Chen et al. [23] conducted accurate simulations of rotating stall propagation for the NASA CC3 compressor with and without the use of an air injection. They used the reverse tangent injection, with a flow rate of 1.7% of the design inlet mass flow rate. Results showed that an air injection can minimize stall areas and prevent reversed flow at the diffuser inlet, and can also correct the flow angle at the impeller exit to be close to the designed value. Lu et al. [47] introduced numerical and experimental research about using micro air injections for controlling stall in an axial compressor. This research pointed to the role of air injections in delaying stall propagation by shifting the tip leakage vortices away from the impeller tip inlet region. The kinetic energy added by the air injection at the tip region was found to be the reason behind enhancing the flow distribution at stall. Khaleghi et al. [48] performed numerical simulations of the air injection parameters in order to clarify the effect of changing the air injection velocity on the compressor stability. Based on the results of this numerical analysis, the maximum gain in the surge margin happens if the Mach number of the injected flow is adjusted to unity. Hiller et al. [24] measured the flow rate corresponding to surge when using air injections with various injection flow rates and temperatures at the first stage of a multistage compressor. The changing of the injected air temperature was found to have little effect on the compressor stability. On the other hand, the injection flow rate of at least 2% of the compressor inlet flow rate improved the compressor’s stable zone. The air injection study of Skoch [22] is valuable in understanding the CC3 compressor, as for the study by Chen et al. [23]; however, both of these studies presented the air injection effect in the forward and tangent directions only, which is inadequate for realizing the optimum injection angle. Furthermore, the injection flow rate and temperature parameters were taken into 16 consideration but without studying the injection angles variation for the study of Hiller et al. [24]. There are few studies which consider all the injection parameters, but the majority of these studies are for axial compressors, and the others do not show a complete optimization of these parameters. In an attempt to find the best injection angle for the CC3 compressor, the variation of the injection flow rate and the injection angle were studied in the presented thesis. 2.3 Structural Analysis of Compressor Blades The compressor blades are affected by aerodynamic fluctuating forces during stall, and these forces cause large deflection and stresses on the blades. For this reason, it is very important to study the effect of the unstable flow condition on blade structure. For any aero elastic system, there are fluid and structure domains affecting each other. The fluid domain causes aerodynamic forces on the structure domain where the stresses and deflections take place. The boundary movement between the two domains depends on both the deformation and the reaction force at the same time. The coupling between the fluid and structure domains is known numerically as the Fluid-Structure Interaction (FSI). The coupling between fluid and structure domains may be weak as in cases of turbo machinery flutter problems, when the blades have a large inertia and the forces acting on it are relatively small [49]. For the case of small blade vibrations, the flow can be simulated by creating a superposition of the steady state flow and a small perturbation flow [50]. The nonlinear flow may result from the small blade vibrations if the operating flow condition is at the transonic region [51, 52]. Gnesin et al. [53] used the time marching method to simulate the blades’ vibration and the aero elastic action for an axial flow turbine. The mode shapes were combined linearly to represent the motion of the blades. Results showed that there is a nonlinear behaviordue to the interaction between the mode shapes. This nonlinearity caused a limitation of the blade vibration oscillations. Srivastava et al. [54] created a numerical model based on solving the Navier Stokes Equations (N.S.E.). The results of the model were compared to the results of other models which deal with a linearized potential flow equation and a linearized Euler equation in order to formulate the blade structural analysis for a turbine and for a fan. The comparison showed that the proposed model as well as the other models predicted the flutter mode and the inter blade 17 phase angle with high accuracy. Also, it was found that the aerodynamic damping decreases as the back pressure increases. Im et al. [55] presented a numerical investigation about the causes of the non-synchronous vibration of axial compressor blades by creating a coupled interaction between the fluid and structure domains. It was found that at the level of 78% of the span at the impeller leading edge, the vorticity and turbulence are very high and there are strong unstable tornado vortices which cause the vibrations of the blades. These tornado vortices are formed due to the complex flow resulting from the interaction between the main flow and the secondary flow, including the tip clearance flow. Brandsen et al. [56] introduced a numerical study about the simulation of the impeller blades’ vibrations for an axial compressor by using an accurate FSI model. The impeller blades were excited with a specific vibration mode and frequency. Results showed that the forces and the blades’ vibration frequencies can be calculated accurately using the transient two-way FSI model. Sadeghi and Liu [57] found that the flutter frequency at low mass flow rate is not the same as the blades’ natural frequency when using the coupled fluid-structure method; this result is different from that of the uncoupled method for which the blades are assumed to oscillate at their natural frequencies. Pei et al. [58] compared the accuracy of the coupled and the uncoupled fluid-structure numerical models to predict the unsteady flow behaviorin a single stage centrifugal pump. Results showed that the coupled model was able to capture the variations of the hydrodynamic force more accurately than the uncoupled model. Results indicated that the pressure values are higher for the coupled model compared to the uncoupled model at any impeller position. At the end of this section, the importance of performing structural analysis for axial and centrifugal compressors is clarified. For the current research, structural models were set with different bearing stiffness values in order to identify the degree of stability at the stall condition in cases where the compressor is controlled with air injection, and also where the compressor is not controlled. 2.4 Synopsis of the Literature Table 2-1 shows the research outcome and the comments and the related proposed research in the current work for the most important studies mentioned in the literature. 18 Table 2-1 Literature summary Title Research outcome Author/ year Comments & proposed research in current work Rotating stall characteristics The tip leakage flow and back flow impingement are the main features of rotating stall in axial compressors. Vo et al. [25] / 2008 It is important to study these two features in centrifugal compressors where the flow nature is more complex. Tip leakage vortex development with time during stall is the main cause of surge appearance. Tomita et al [30] / 2013 The tip leakage vortex motion will be reported in order to introduce better understanding of the surge initiation. Tip leakage flow results from the formation of vortices with tornado shape due to flow separation. Yamada et al [29] / 2013 The flow separation locations will be studied in more details. Casing treatment Holed casing treatment can increase the compressor surge margin by enhancing the performance at low and high flow rates. Xu et al. [20] / 2011 As a comment, this is a successful method but there are manufacturing difficulties for producing these holes. The most successful location to place the casing grooves in axial compressors is at the mid of the chord distance. Ivor Day [39] / 2012 This result will be verified for the case of centrifugal compressor by testing three different locations for the grooves. The casing grooves performance can be increased by changing their dimensions. Kim et al. [37] / 2013 This analysis will be expanded in the current research by changing the aspect ratio of the grooves. Air injection The surge margin of a high speed centrifugal compressor was increased by using forward and reverse air injection. Skoch [22] / 2003 The effect of using various air injection angles on the compressor stability will be studied. The maximum gain in the compressor surge margin can be obtained if the injection air Mach number equals one. Khaleghi et al. [48] / 2008 As a comment, this is a very important result which clarifies the optimum injection velocity. Numerical simulations proved the effectiveness of using reverse tangent air injection in improving the compressor stability. Chen et al. [23] / 2009 As a comment, this study was beneficial for knowing the numerical capabilities in simulating air injection. Structural analysis The coupling between fluid and structure domains may be weak as in cases of turbo machinery flutter problems. Sadeghi, M., and Liu, F. [49]/ 2005 As a comment, the large inertia of the blades supports this point and the d blades deflections are relatively small. The causes of the non-synchronous vibration of axial compressor blades were investigated by creating a coupled interaction between the fluid and structure domains. Im et al. [55] / 2012 The vibration and deformation of blades will be studied in the case of stall/surge. The accuracy of the coupled and the uncoupled fluid-structure numerical models to predict the unsteady flow behavior in a single stage centrifugal pump was illustrated. Pei et al. [58] / 2013 This study clarified the difference between coupled and uncoupled fluid-structure interaction simulations and presented a guidance to identify the accuracy for each method. 19 3)CHAPTER 3: NUMERICAL METHODOLOGY This chapter discusses the numerical methodology and the governing equations for the modeling of the turbulent flow. In this research, the commercial program FLUENT [59] was used to perform the CFD simulations using the finite volume method. In the finite volume method, the partial differential equations of continuity, turbulence, energy, and Navier Stokes equations are discretized into a set of algebraic equations. Then, the algebraic equations are linearized and solved to obtain the values of the dependent variables. The 2nd order discretization scheme was applied for space and time for the purpose of increasing the solution accuracy. The conservation equation of transport and its discretization are shown in the Appendix (A.1). The governing equations were solved using the density based method [60]. The solution procedure for the density based method begins by obtaining velocity components from the equations of momentum; then, the velocities are used to compute the density from the equation of continuity, followed by obtaining the pressure terms from the equation of state. The final step is to solve the turbulence equations as indicated in Fig. 3-1. The convergence is achieved for each time step by repeating the steps of this solution procedure. Figure 3-1 Overview of the density based method. 20 3.1 Turbulence Modeling There are a variety of choices for modeling turbulent flow as ݇ െ ߝ , ݇ െ ߱, Reynolds stress, and large eddy simulation models. The main factors which affect the choice of the suitable turbulence model are the physics of the problem, the required accuracy level, computational resources and total simulation time. The realizable ݇ െ ߝ model was chosen to simulate the turbulent flow because it provides results with relatively high accuracy when simulating rotating flow at high speeds and also boundary layer and complex flows [61, 62]. Another reason for selecting the realizable ݇ െ ߝ model is that it does not have the major disadvantages of the other turbulence two equation models such as the standard ݇ െ ߝ model which provides inaccurate results for complex flow with swirl, and the RNG ݇ െ ߝ model, which is subjected to limitations due to the isotropic eddy viscosity assumption. Also, the results of the realizable ݇ െ ߝ model have an acceptable accuracy and can be achieved with less computational time compared to the Reynolds Stress Model. The mean velocity gradients are related to the Reynolds stresses according to the Boussinesq hypothesis [63] as indicated in the Appendix (A.2). The exact transport equation of the mean square vorticity fluctuation in the realizable ݇ െ ߝ model is used to derive the turbulence dissipation rate. The turbulent viscosity is a function of the velocity gradients and is not constant. Many comparisons were performed between various turbulence models to verify the accuracy of simulating the distribution of pressure and velocity. These comparisons proved the high prediction accuracy of the realizable ݇ െ ߝ model for unsteady high speed flow. Details about the realizable ݇ െ ߝ equations are presented in the Appendix (A.3). By using the realizable ݇ െ ߝ turbulence model, the turbulent viscosity prediction is enhanced compared to other models as the standard ݇ െ ߝ model by taking advantage of using the term ܥఓ as a function of velocities and shear stresses instead of being constant. 3.2 Near-Wall Treatments for Turbulent Flow At the areas very close to the wall, the tangential velocity fluctuations are reduced by the effect of viscous damping, and the normal velocity fluctuations are affected by kinematic blocking. The near-wall area is composed of three regions; the viscous sub layer, the interim region and the fully turbulent layer. Figure 3-2 shows the layers of the near wall region plotted in 21 semi-log coordinates. The vertical axis in Fig. 3-2 represents the dimensionless velocity (ܷା) while the horizontal axis is for the dimensionless wall distance (ݕା). Definitions of the near-wall treatment parameters are mentioned in the Appendix (A.4). Figure 3-2 Classification of the near-wall region. For the viscous sub layer, the flow is nearly laminar and is governed by molecular viscosity which has the greatest effect on momentum and heat transfer. At the outer layer or the fully turbulent layer, the turbulence plays a dominant role. The mean velocity has large gradients close to the outer layer and this increases the turbulence’s kinetic energy and turbulence rate. The turbulence and the molecular viscosity have nearly the same impact on the flow at the region between the viscous sub layer and the fully turbulent layer [64]. 3.2.1 Wall Functions Wall functions are a combination of semi-empirical formulas which work as a link between the flow variables at the wall and the corresponding values at the near wall cells. The 22 idea of the standard wall functions was firstly discussed by Launder and Spalding [65]. There are wall functions for momentum, energy and turbulence. (a) Using wall functions for momentum The general concept of the near wall treatment for the turbulent flow was discussed before in the introduction of section 3-2 based on the definitions of ݕାand ܷା. For most of the CFD programs, the same concept is applied but using other definitions. Instead of using ݕାand ܷା, the dimensionless quantities ݕ∗ and ܷ∗ are used as dimensionless distance and dimensionless velocity, respectively. The quantities ݕ∗ and ݕା are nearly equivalent for the equilibrium turbulent boundary layers. The law for the mean velocity is valid for (300 >ݕ∗>30). The log-law is used in the CFD program [59] when ݕ∗> 11.225. More details about the equations for the using of wall functions for momentum can be found in Appendix (A.5). (b) Using wall functions for energy As for the momentum calculations where the logarithmic law for the mean velocity is used, there is also another logarithmic law applied for mean temperature. There are two different wall laws for the temperature; the linear law and the logarithmic law. The linear law is applied when the conduction is important and the logarithmic law is activated when the turbulence effects are dominant. The thermal conduction layer thickness (ݕ∗் ) is not the same as the viscous sub layer thickness and also differs from one fluid to another. In the case of compressible flow at high speeds, the temperature distribution close to the wall is different from that for the subsonic flows due to the heating caused by viscous dissipation [66]. The law for temperature and the procedure of solution are described in the Appendix (A.6). (c) Using wall functions for turbulence For the modeling of turbulence using the ݇ െ ߝ model, the turbulence kinetic energy equation is solved in all domains including the near wall areas. Equation 3-1 shows the boundary condition applied at the wall for the turbulence kinetic energy. (3-1) Where n is a local coordinate normal to the wall. 23 Based on the local equilibrium hypothesis, the kinetic energy production term ܩ and its dissipation rate ߝ are calculated at the near wall area as described in equations 3-2 and 3-3, respectively. The turbulence dissipation rate equation of the ݇ െ ߝ model is not solved but it is calculated from equation 3-3 instead. (3-2) (3-3) Finally, it may be seen that the wall boundary conditions for the solution variables including ݇, ߝ, velocity and temperature are all considered by the wall functions. 3.3 Dynamic Mesh Model The dynamic mesh model is used for modeling flows which contain moving and static domains and it may also be used to simulate the flow when the domain shape changes with time due to the motion of the boundaries. 3.3.1 Dynamic Mesh Concept The overall mesh contains rotating domains and fixed domains and it is important to simulate the interaction between the rotating and fixed domains accurately especially at high speeds. For this reason, the dynamic mesh model was used and its idea is to make the rotating domain slides on the fixed domains at the interface surfaces which separate these domains. At the interface, the mesh must be modified and reformed at the end of every time step. At each time step, the intersection between the interface zones is calculated and then the interface flux is determined. For periodic problems, this intersection results in an interior zone if the two interface zones overlap, and also results in periodic zones if there is no overlap. The interior zone has fluid cells on both sides of the intersection. For non-periodic problems, wall zones are produced instead of the periodic zones. Figure 3-3 shows an example for a non-periodic flow problem. Two objects are moving in opposite directions in a fluid medium (Fig. 3-3 (a)). In this example, the whole fluid domain is split into two zones; each zone covers one 24 moving object and slides over the other zone through the interface zones (Fig. 3-3 (b)). During the sliding process between the two zones, the common area is called the interior zone, which is surrounded by fluid from both sides, and the remaining areas are treated as wall surfaces (Fig. 3-3 (c)). Figure 3-4 shows an example of the faces created due to the interaction between interface zones. The interface zones are composed of four different faces; A-B, B-C, D-E and E-F. The faces a-d, d-b, b-e, e-c and c-f are the result of intersection between the interface zones. The overlapping area between cell zone 1 and cell zone 2 represents the interior zone which includes faces d-b, b-e, and e-c. The other faces outside of the overlap area, including faces a-d and c-f, represent the periodic zone. The flux is computed across the interface into the cell III, for example, by ignoring face B-C and using faces b-e and e-c instead by transferring the information from cells IV and VI into cell III. (a) Two objects moving in opposite directions (b) Splitting the domain with the formation of two interfaces (c) Description of the interior and wall zones formed during objects motion Figure 3-3 Schematic drawing for the various zones result due to a non-periodic interaction. 25 Figure 3-4 Faces and zones created by the interaction between two cell zones. 3.3.2 Dynamic Mesh Conservation Equations The following equations explain the main concept of the dynamic mesh theory [60, 67,68]. For a general flow variable ∅, the governing equation used for the sliding mesh method is described in equation 3-4. ݀݀ݐ න ߩ∅ܸ݀ න ߩ∅ሺݑ െ ݑడሻ. ݀ܣ ൌ න Γ∅. dAడ න S∅dV (3-4) Where V is the volume of the mesh and ߩ is the density. While ݑ and ݑ are the fluid velocity vector and the velocity of the moving grid mesh. Also, S∅ is the source term of ∅ and Γ is the diffusion coefficient. For second order accuracy, the volume integration with time is calculated as shown in equation 3-5. ݀݀ݐ න ߩ∅ܸ݀ൌ 3ሺߩ∅ܸሻାଵ െ 4ሺߩ∅ܸሻ ሺߩ∅ܸሻିଵ2∆ݐ (3-5) 26 Where n represents the current time step and the term ܸାଵ is calculated as described in equation 3-6 by relating it to the volume at the current time step and the rate of change of volume with time. ܸାଵ ൌ ܸ ܸ݀݀ݐ ∆ݐ (3-6) The rate of volume change with time (ௗௗ௧) is computed by integrating the grid velocity over the area (equation 3-7). ܸ݀݀ݐ ൌ න ݑడ. ݀ܣ (3-7) As an equivalent method to the integration, the term ௗௗ௧ is computed by making a summation of the dot product of the velocity vector and area vector at each face of the mesh cell volume as shown in equation 3-8. ܸ݀݀ݐ ൌ ݑ. ܣ (3-8) Where ݊ represents the number of faces for each control volume and ܣ is the area vector for the face j. The R.H.S. of equation 3-8 is a summation of the volume change with time for each cell or in other words, ݑ. ܣ ൌ ఋೕఋ௧ which represents the volume swept by the control volume j with time. The value of the term ݑ. ܣ at a specific time step is determined by computing the same term at previous time steps according to equation 3-9. ሺݑ. ܣሻାଵ ൌ32 ሺݑ. ܣሻ െ 12 ሺݑ. ܣሻିଵ ൌ 32 ሺߜ ܸߜݐ ሻ െ 12 ሺߜ ܸߜݐ ሻିଵ (3-9) 3.4 Fluid-Structure Interface Fluid structure interaction (FSI) is a definition for problems which have a strong dependency between solids and fluids. There is what is called an interface between the solid and the fluid domain. At the interface, all governing equations and also the boundary conditions for the fluid and solid regions must be satisfied simultaneously. For the purpose of connecting the governing equations from the fluid and solid domains, a set of coupling conditions are applied at 27 the interface. There are two types of coupling conditions; the kinematic and the traction (dynamic) conditions. The kinematic condition is related to the motion where the fluid and the structure must have the same velocity at the interface surface while the traction condition satisfies the force balance and insures that the fluid and solid are in local equilibrium. Due to the difference in the accuracy requirements for the fluid and the structural simulations, there are different surface discretizations used. Usually, the fluid domain requires a finer mesh comparing with the solid domain in order to capture the nonlinear features accurately. This discretization difference causes mismatch of the mesh points at the interface surface. Due to the existence of the mesh mismatch, there is a numerical approach to transfer the pressure values on the mesh of the fluid domain to corresponding nodal forces on the solid domain and also for the transfer of displacements of nodes at the solid domain to the corresponding perturbations on the fluid domain. In order to ensure the consistency and the conservativeness of the load transfer and the displacements transfer processes, the conservation of load and energy must be satisfied where the sum of the nodal forces in the solid domain must be equal to the integral of the pressure fields on the surface boundary in the fluid domain. Regarding the numerical discretization, there are two approaches; the monolithic (fully coupled) approach and the partitioned (staggered) approach. For the monolithic approach, the governing equations of the fluid and solid are solved together and this requires a conformal mesh with matching nodal positions between the fluid mesh and the solid mesh and also requires fully integrated FSI solver. For the partitioned approach where the mesh can be non-conformal, separate solvers can be used to solve the fluid and the solid phases, in sequential order. Also, for the partitioned approach, there are two solution methods namely one-way coupling and two-way coupling. The one-way coupling was used in the current research because of the relatively small deformation of the blades. In the one-way coupling method, the fluid pressure on the fluid domain is transferred to the solid domain, but the displacements are not transferred back to the fluid domain. For the data transfer of System Coupling in the current research, when conservative quantities as mass and momentum are transferred, a method called General Grid Interface is applied. This method calculates weights that make the transferred quantity conserved locally. If all nodes are mapped correctly, the method is also globally conservative. 28 4)CHAPTER 4: SELECTED CASE STUDY MODELING AND VALIDATION 4.1 Description of the Case Study The NASA CC3 high speed centrifugal compressor was selected for this research on rotating stall simulations because it has been involved in many previous experimental and numerical studies [22, 23, 46,69-72]. Figure 4-1 (a) and Fig. 4-1 (b) show the cross section view and the three dimensional view of the NASA CC3 compressor, respectively. It can be seen in Fig. 4-1 that the compressor impeller is composed of 15 main blades and 15 splitter blades, while the diffuser contains 24 vane island shaped passages. Table 4-1 summarizes the main dimensions of the impeller and diffuser and the specification of the compressor design condition. At its best efficiency point, the design conditions are 21,789 rpm rotational speed, and a mass flow rate of 4.54 kg/s, with a total to total pressure ratio of 4. McKain and Holbrook [73] reported the complete design specifications of the NASA CC3 compressor, along with a vibration and structural analysis of the compressor blades. (a) Cross section view (b) 3-D view Figure 4-1 NASA CC3 centrifugal compressor [74] Impeller Diffuser29 Table 4-1 NASA CC3 compressor specifications Compressor geometry Impeller Inlet tip diameter 210 mm Exit tip diameter 431 mm Number of main blades 15 Number of splitter blades 15 Tip clearance at impeller exit 0.2 mm Diffuser Shape Vane Island Inlet diameter 465.5 mm Exit diameter 724 mm Number of vanes 24 Design operating condition Pressure ratio (total to total) 4 Adiabatic efficiency (total to total) 83.30% Rotational speed 21,789 rpm Mass flow rate 4.54 kg/s 4.2 Numerical Modeling of the Case Study 4.2.1 CFD Models There are two numerical models used to simulate the flow in the NASA CC3 compressor: the first is the single passage model (Fig. 4-2) and the second is the full passage model (Fig. 4-3). The single passage model was used only for the purposes of validating the numerical results with measurements while the full passage model was used in all other simulations as rotating stall and its control. The idea of the single passage model is to simulate part of the compressor volume to save computational time, and it is composed of three parts: the inlet guide duct, the impeller, and the diffuser, including the exit bend as shown in Fig. 4-2. This model was used in the steady state simulations where the flow distribution does not change from one passage to another. The full annulus model contains the same three main parts as the single passage model does, but also includes the whole compressor volume as indicated in Fig. 4-3. The full annulus model was used in the unsteady simulations in order to capture the flow variations accurately. 30 Figure 4-2 The single passage numerical model Figure 4-3 The full annulus numerical model There are many boundary condition types available for the numerical simulations at the flow inlet and outlet, like velocity inlet, mass flow rate, pressure inlet, outlet flow and pressure exit boundary conditions. Some of these boundary conditions are not suitable for high speed Bend Full impeller blade Splitter impeller blade DiffuserX-directionY-directionAxial direction 31 compressible flow. For example, in the velocity inlet boundary condition, the density and velocity have fixed values at the inlet and this causes error in the compressible flow simulation. Another example is the mass flow rate boundary condition, which does not provide the required accurate pressure distribution. For the purposes of obtaining a fixed pressure ratio between inlet and outlet, the pressure inlet boundary condition was assumed at the compressor inlet and the pressure outlet boundary condition was used at the compressor exit. In this way, by fixing the pressure values at inlet and outlet, the solution converges to the corresponding flow rate. The total pressure of 101,325 pa and total temperature of 288 K with a turbulence intensity of 10% were specified as inlet conditions according to the available data and measurements [73]. The static pressure was specified at the exit boundary with a value depending on the simulated operating condition. For the air injection simulations, the mass flow rate boundary condition was used at the injection inlet locations in order to have the ability to enter the injection flow rate with a specific injection angle. The sliding mesh was used, as explained in Chapter 3, to simulate the interaction between the rotating zones, including the impeller, and the stationary zones, including the diffuser, through interface surfaces. 4.2.2 Structural Model The numerical model used for the structural analysis represents the impeller solid volume (Fig. 4-4). The global algorithm mapping method was implemented in order to transfer stresses from the fluid domain to the structural domain [75]. Figure 4-4 The structural model 32 At some locations, there is a gap between the two meshes; in this case, the node on the fluid mesh is mapped to the closest element on the structural mesh. Figure 4-5 (a) shows that for the node N1, there are two facing elements e1 and e2; this node is mapped to the element e1 because the distance d1 is less than the distance d2, and this means that the selection criteria is based on minimizing the mapping transfer distance in the gap between the two meshes. In the case that there is any misalignment between the interface edges, the mapping is done by selecting the closest node as shown in Fig. 4-5 (b). (a) Node mapping according to the closest distance to the facing elements (b) Node is mapped to the closest facing node for the case of misaligned edges Figure 4-5 Global mapping method [75] The zero displacement boundary condition was applied at a specific location on the shaft near the inlet. The same boundary condition was also applied at another location on the shaft near the exit; however, at this location, the displacement in the axial direction is left free to simulate a sliding joint. The above boundary conditions matched those used for the NASA model and, therefore, the results were used for verification and validation of the current model [73]. 4.2.3 Meshing Analysis The mesh sensitivity analysis was performed by comparing the results of the numerical model when using various mesh sizes. The total number of mesh elements in this comparison ranges from 5 million to 15 million cells. This variety in the total number of mesh elements was achieved by changing the size of the smallest mesh element and the growth rate of the mesh at different locations. The following findings present the comparison between the results of three mesh set levels; the first level represents the coarse mesh with 5 million cells, while the second and third levels represent mesh with 10 and 15 million cells, respectively. Figure 4-6 (a) shows 33 the pressure variation along the axial direction at the mid of the span length for the three mesh levels with the calculation of the extrapolation curve which was obtained using the method of Richardson [76]. Figure 4-6 (b) indicates the variation of the pressure in the case of mesh with 15 million cells, and presents the error bars which represent the GCI values. The error values range between 0.65% and 2.3%. These error values are relatively low and indicate that the 15 million cell level is sufficient for performing the simulations with high accuracy. (a) For three different mesh cases and for the calculated extrapolation curve. (b) For the fine mesh case clarified with error bars. Figure 4-6 Pressure variation at the mid span with the axial direction (z) The ratio between the time step and the size of the grid was altered in order to test the temporal accuracy. For the unsteady simulations, the time step used was 6x10-6 second because it was found that the solution convergence can be achieved within only 20 iterations per time step, by using this time step value. Also, the CFL number was varied for the selected time step in order to study its effect on the solution stability. Figure 4-7 indicates the mass flow rate variation with the number of iterations for three different values of the CFL number. The CFL number and time step have a direct effect on the solution convergence. It may be concluded from Fig. 4-7 that the smallest fluctuations and the fastest convergence rate happen for the case where CFL 34 equals 3. The convergence is reached after nearly 30,000 iterations at the compressor design operating condition corresponding to a rotational speed of 21,789 rpm and a pressure ratio of 4. Figure 4-7 Mass flow rate convergence for different CFL values The air injection method was used in this research to control stall; in this method, the air is injected at certain points in the vaneless region as explained in the results chapter. The mesh at the vaneless region was constructed by taking into consideration the mesh quality at the injection locations. Each injection hole has a diameter of 2 mm and the mesh at this hole has a starting size of 0.2 mm and increases gradually with a growth rate of 1.1 (Fig. 4-8). Figure 4-8 Meshing at the vaneless region 30◦ Air injection hole diameter = 2 mm Starting mesh size of 0.2 mm with mesh growth rate of 1.1 R = 0.234 m 35 The mesh was constructed in a way that it is fine enough at the areas of small gaps, as in the tip clearance region, in order to obtain accurate simulation of the flow movement above the blade tips. The tip clearance value is 0.2 mm, and the starting mesh size used is 0.05 mm, representing the thickness of the first mesh layer at the shroud surface. A mesh growth rate of 1.1 was applied in order to increase the mesh size as it moves away from the shroud surface. The mesh was also refined at the blade surface and at the leading and trailing edges of the impeller blades for better prediction of velocity gradients at these surfaces. Figure 4-9 shows the final mesh generated for the CFD numerical model with a total number of mesh elements of 15 million cells. On the other hand, for the meshing of the mechanical numerical model, several trials were performed to optimize the mesh quality by monitoring the mapped variables from the CFD mesh. The final mechanical model mesh is about 2 million cells as indicated in Fig. 4-10. Figure 4-9 Meshing of the CFD model Figure 4-10 Meshing of the mechanical model 36 4.3 Validation of the Numerical Model Different operating conditions were simulated in order to form the pressure and efficiency maps of the compressor at 70%, 80%, 90% and 100% of the compressor design speed. Figure 4-11 (a) and Fig. 4-11 (b) show the pressure map and the efficiency map respectively for the numerical model results plotted with measurements done by Skoch [69]. The numerical model provided correct predictions of the pressure and efficiency curve trends at different speeds. There is a relatively small error between measurements and numerical results in the range of 2% to 5%. All simulations were done using a high speed computer with parallel processing (16 cores, 2.1 GHZ, 96 GB Memory). For the steady state simulations using the single passage model, the complete solution converged after 30,000 to 40,000 iterations using the time step of 6x10-6 second, with a total computational CPU time of nearly 40 hours for each simulation case. (a) Pressure map (b) Efficiency map Figure 4-11 Compressor pressure and efficiency maps (numerical results versus measurements [69]) The numerical accuracy was also verified by comparing the pressure and Mach number contours at blades surfaces at the design operating condition of the presented numerical model 37 with those from a previous numerical study [75] as indicated in Fig. 4-12. The pressure distribution of the current numerical model (Fig. 4-12 (a)) is close to the pressure distribution found in the previous study (Fig. 4-12 (b)) because the values and locations of the minimum and maximum pressure regions are similar for both models. The Mach number contours from the previous study (Fig. 4-12 (d)) have slightly higher values at the inlet tip compared to the present numerical model results (Fig. 4-12 (c)) but the main variations along the span length are the same for the two models. (a) Pressure contours [for the current numerical study] (b) Pressure contours [for a previous numerical study] (c) Mach contours [for the current numerical study] (d) Mach contours [for a previous numerical study] Figure 4-12 Pressure and Mach number contours for the current model results and for previous CFD model results [75] Figure 4-13 indicates the average velocity variation with the percentage meridional length for the numerical model results and also for previous measurements [73]. It may be concluded from Fig. 4-13 that there is a good match between the two curves in terms of the 38 slope, except for a small margin between 30% and 40% of the meridional distance where the numerical results curve has lower values. Figure 4-13 Mean velocity variation with the % meridional distance (numerical results versus measurements [73]) It is also important to check the accuracy of velocity distribution from hub to tip. The most complicated flow variations occur at the vaneless region especially near the shroud surface. Figure 4-14 shows the distribution of the velocity relative to the tip value from hub to tip at a point inside the vaneless region at a radius ratio of 1.1 for the numerical model results and for measurements [69]. Figure 4-14 Velocity profile from hub to tip at the vaneless region at radius ratio of 1.1 (numerical results versus measurements [69]) 39 It may be noted from Fig. 4-14 that the numerical model successfully predicted the velocity profile along the span length with high accuracy. Figure 4-15 shows the variation of the static pressure normalized by the total inlet pressure, with the normalized meridional chord at the impeller tip for the blade suction and pressure surfaces. The pressure variation shown in Fig. 4-15 is for the numerical model results compared with measurements [77]. It may be concluded from Fig. 4-15 that the blade pressure side curve predicted by the numerical model matches the measurements accurately with a small deviation at the leading and trailing edge points only. There is also a relatively small error between the numerical results and measurements at the mid-range of the suction side curve, but in general, the blade loading curves were predicted correctly with a maximum relative error of 5 %. Figure 4-15 Normalized pressure variation with the normalized meridional distance at the impeller tip (numerical results versus measurements [77]) Validation for the flow in the unsteady case was performed in order to ensure the accuracy of the CFD code for predicting transient flow variations. Figures 4-16 (a) and 4-16 (b) show the surge cycle for measurements [72] and for the numerical model results of this study, respectively. When comparing one of the measured surge cycles shown in Fig.4-16 (a) to that in Fig.4-16 (b), there are some fluctuations in the beginning with nearly constant amplitude and at some time the pressure turns up to approximately 400 kpa and then decreases sharply. 40 (a) Surge cycle curve [measurements] (b) Surge cycle [numerical results] (c) Point A [measurements] (d) Point A [CFD results] (e) Point B [measurements] (f) Point B [CFD results] (g) Point C [measurements] (h) Point C [CFD results] Figure 4-16 Velocity at the diffuser different times of the surge cycle (numerical results versus measurements [72] Point APoint B Point C41 The numerical results are close to the measurements especially for the surge cycle time interval and the fluctuations amplitude. Figure 4-16 ((c) to (h)) shows the velocity distribution enhanced with velocity vectors inside the diffuser at four different points, as clarified on the surge cycle diagram for the numerical results of this study versus experimental results of reference [72]. At the steady state condition (Point A), the comparison between Fig.4-16 (c) and Fig.4-16 (d) indicates that the flow enters the diffuser with a maximum velocity value around 375 m/s, and the flow distribution for the CFD results is consistent with the measurements. At the maximum pressure point of the surge cycle (Point B), the CFD results (Fig.4-16 (f)) confirm the measurements (Fig.4-16 (e)), revealing that there is a reversed flow inside the diffuser with a maximum velocity of 225 m/s at the diffuser throat. Finally, at the minimum pressure value (Point C), the measurements (Figure 4-16 (g)) show that a small portion of the flow domain with a corrected flow direction appears at the diffuser throat, while the CFD results (Figure 4-16 (h)) confirm this finding but with a wider corrected flow region. It may be concluded from Figure 4-16 that CFD modeling can capture the unsteady flow variations during surge accurately. One of the most important parameters needed for the validation of the mechanical model is the stress distribution. Figure 4-17 (a to c) shows the maximum principal stress at the full blade pressure surface, the equivalent stress at the full blade pressure surface, and the equivalent stress at the splitter blade pressure surface, for the results of the currently presented numerical model. Figure 4-17 (d to f) indicates the same parameters as Figure4-17 (a to c), but for the numerical model results done by NASA [73]. By comparing Figure 4-17 (d) with Figure 4-17 (a), the currently presented numerical model results agree with the NASA model results that there are two regions with high stresses within the range from 27 KSI to 31 KSI at the hub and shroud surfaces close to the impeller exit, and the location with the minimum amount of stress is at the blade inlet. The equivalent stress distribution shown in Fig. 4-17 (e) demonstrates that the region of low stress at the blade tip inlet is expanded for a large area and the maximum stress appears at two areas with a value of 30 KSI; these results are very close to those shown in Fig. 4-17 (b). Figure 4-17 (f) indicates that the maximum stress regions are enlarged compared with those in Fig. 4-17 (e) and the minimum stress region is localized to a relatively smaller region; this result is confirmed by the NASA model results shown in Fig. 4-17 (c). 42 (a) Maximum principal stress at the full blade pressure surface (for NASA results [73]) (b) Equivalent stress at the full blade pressure surface (for NASA results [73]) (c) Equivalent stress at the splitter blade pressure surface (for NASA results [73]) (d) Maximum principal stress at the full blade pressure surface (for the presented numerical model results) (e) Equivalent stress at the full blade pressure surface (for the presented numerical model results) (f) Equivalent stress at the splitter blade pressure surface (for the presented numerical model results) Figure 4-17 Stress distribution at impeller blade surface (The presented numerical model results compared with the NASA numerical model results [73]) The deflection of the blades was computed in the axial direction as presented in Fig. 4-18. The most critical point is at the blade tip exit, where the axial deflection is - 0.014 inch according to the numerical model results, while the corresponding value for the NASA model results is - 0.01463 inch. This point at the blade tip exit is critical because it is the closest point to the shroud surface compared to any other point on the blade tip surface [73]. This result may be considered as a validation for the deflection calculations. PSI PSI PSI 43 Figure 4-18 Axial deflection (for the presented numerical model results) The most critical point (at the blade tip exit) in terms of its closeness to the casing +ve axial direction44 5)CHAPTER 5: RESULTS AND DISCUSSION 5.1 Rotating Stall Simulation for the Compressor with a Vaned Diffuser 5.1.1 Overall View of Stall Characteristics The development of stall in the early stage of surge may be viewed in the form of main characteristics of the stall phenomenon. Due to weakness of the inlet axial flow velocity, separation occurs early, causing turbulence and vortices in the impeller passages as shown in Fig. 5-1, in which the velocity vectors are presented at a plane close to the impeller blades tip during the stall process. This early separation can be seen in the magnified view in Fig. 5-1, where the flow separation appears at the impeller’s leading edge, followed by a stalled region. Early separation causing vortex formation Figure 5-1 Velocity vectors near impeller shroud during stall The first main characteristic of the stall phenomenon is the occurrence of the tip leakage flow at which the fluid moves from blade pressure side to blade suction side across the clearance between the blade tip and the shroud, causing disturbance in the flow in the affected passage (Fig. 5-2). 45 Figure 5-2 Velocity vectors showing tip leakage flow The second main characteristic of the stall phenomenon is the back flow impingement when the fluid at the rotor exit is turning back from one passage to another through the vaneless region instead of transferring into the diffuser, and this causes fluid velocity deviation at the impeller exit as shown in Fig. 5-3. Figure 5-3 Velocity vectors showing back flow impingement at impeller exit during stall The velocity deviation at the impeller exit causes diffuser blockage or failure of the fluid to enter the diffuser at the designed vane angle; and this causes the formation of very low velocity regions at the diffuser entrance, which can develop over time to form a back flow: this is considered the sign of the beginning of a surge. Figure 5-4 shows the velocity contours near the 46 diffuser shroud surface during stall, and it is clear that there are some diffuser passages with very low velocity values. It is also shown in Fig. 5-4 the velocity vectors at part of the compressor (the area inside the rectangle); and the velocity vectors indicate the formation of back flow, which explains the low velocity values at this area. Velocity vectors show formation of back flow Figure 5-4 Velocity contours near diffuser shroud surface during stall 5.1.2 Tip Leakage Flow Variation With Time The path lines are the trajectories of flow particles starting from the inlet through the impeller and diffuser. Figures 5-5 (a) to 5-5 (d) show path lines colored by velocity magnitude at different times (2, 4, 7 revolutions)at a condition near surge (reference time = 0). The path lines are slightly disturbed from the ideal path and there are no odd flow patterns in the vaneless region or in the diffuser (Fig. 5-5 (a)). 47 m/s (a) At the beginning of stall (reference time=0) (b) After 2 revolutions (c) After 4 revolutions (d) After 7 revolutions Figure 5-5 Path lines near shroud surface at different time frames Appearance of vortex flow Flow direction deviation causes diffuser blockageTip leakage flow increases Appearance of tip leakage 48 Figure 5-5 (b) indicates that there are some path lines crossing the impeller blades through the tip clearance and concentrated near the impeller exit, indicating the initiation of the tip leakage flow. Figure 5-5 (c) shows that the number of path lines crossing the tip is increasing and starting to appear at lower radii, and are not only appearing near the impeller exit. The increase in the tip leakage flow is developed because the impeller passages are partially blocked at the impeller exit. The worst case scenario is the occurrence of reversed flow (Fig. 5-5 (d)), where the flow in the vaneless region is very weak and the flow angle at the diffuser entrance is very high, so the fluid is turned back towards the impeller blades; this is considered to be the beginning of the surge cycle. Figure 5-5 (d) also shows the formation of vortices in some passages close to the vaneless region. Figure 5-6 describes velocity vectors at five planes perpendicular to flow direction in the impeller at times corresponding to 4 and 7 revolutions from the reference time stated previously at an operating condition near to surge line. m/s (a) After 4 revolutions from the reference time (b) After 7 revolutions from the reference time Figure 5-6 Velocity vectors at different planes through the impeller at an instant of time during stall The tip leakage vortex occurs because of the tip leakage flow interaction with the incoming radial inward flow as seen in Fig. 5-6 (a) at some passages near the exit. Also, the tip leakage vortex develops at the splitters inlet as shown in Fig. 5-6 (b). As the tip leakage flow 49 increases and the axial flow velocity component decreases, the tip leakage vortex starts to appear earlier in the impeller near inlet, while the tip leakage vortex near the exit becomes larger. It can be concluded that there are two vortices, one is perpendicular to the flow direction (the tip leakage vortex), and the other is at a plane parallel and close to the impeller shroud. 5.1.3 Rotating Stall Effects on the Compressor Blades During the stall process, the velocity and pressure distributions vary with time, especially near the shroud surface where the stall cells initiated. Figure 5-7 indicates the pressure distribution at a plane of 95% of the span length at different time frames, starting from the reference time at which stall is initiated and after 3, 5, 7 and 9 rotor revolutions from the reference time. It may be noted that pressure distribution is normal at the reference time with a pressure ratio of 4:1, but after 3 revolutions, pressure distribution changes at the vaneless region between the impeller and the diffuser and this can be related to the velocity disturbances. After 5 revolutions, some localized high pressure regions appear at the vaneless region; this is because the fluid is impacting the impeller blades. These high pressure regions increase in size and number, which can be seen after 7 and 9 revolutions. At the reference time After 3 revolutions After 5 revolutions After 7 revolutions After 9 revolutions Figure 5-7 Pressure development during stall 50 When the stall starts, there are pressure and velocity fluctuations at some places in the compressor, especially near the shroud surface; these fluctuations continue and may increase with time to cause surge. Figure 5-8 shows the fluid stream lines at 95% of the span length through the impeller and diffuser, starting from the stall inception condition up to the condition after surge initiation. It can be noted from Fig. 5-8 that the stream lines are normally distributed with little disturbances at the diffuser exit and the vaneless region when the stall initiates. After 5 rotor revolutions from the beginning of stall, flow separations start to appear at the impeller blades’ leading edge and the fluid deviates from the ideal flow path. After 7 rotor revolutions from stall initiation, separations are seen in most of the impeller passages. Furthermore, the stream lines in the vaneless region deviate toward the tangential direction, and there is back flow in some of the diffuser passages. The beginning of surge is considered to occur after 9 rotor revolutions because the number of diffuser passages with back flow increases, and there is a complete breakdown of flow distribution with complete back flow within 2 rotor revolutions from the beginning of surge. At the reference time (stall initiation) After 5 revolutions After 7 revolutions After 9 revolutions (start of surge) After 11 revolutions Figure 5-8 Stream lines distribution during stall and surge at location near shroud surface 51 Pressure fluctuations during stall can be expressed in the form of forces. Figure 5-9 shows the variation of the force acting on the impeller blades in the axial direction with time during stall and surge. Figure 5-9 Variation of the force affecting impeller blades with time during surge There are three curves indicating: the force on the tip of the blades, the force on blades suction and pressure surfaces, and finally, the total force which equals the summation of the forces of the other two curves. It can be concluded from Fig. 5-9 that the total force is increasing up to double of its initial value. Meanwhile the force acting on the blade tips increases by 1.5 times its initial value during the surge event, and this increase in force is followed by another decrease as a part of the surge cycle which can be repeated with time; this means that the blades may be damaged due to this strong variation in force. 5.1.4 The Most Affected Region by Rotating Stall It is important to monitor the flow at the tip clearance, because this provides an indication of the strength of the tip leakage flow. Five points located at the tip clearance gap above the impeller blade tips were selected in order to monitor the dynamic pressure as shown in Fig. 5-10. Figure 5-11 illustrates the dynamic pressure variation with the number of time steps during stall 52 at the five points selected and listed in Fig. 5-10, where the time step used in the simulations is 2x10-6 second. It can be noted from Fig. 5-11 that for all pressure signals, there are fluctuations with time, and that the amplitude of fluctuations varies from one point to another because there is a different pressure distribution for each point according to its location. The fluctuations at all points at a specific time point begin to increase suddenly and this is an indication of stall inception. Signals at points 1 and 2 have very large variations compared to other points. This means that the tip leakage flow at points 1 and 2 or at the impeller blade tips’ trailing edge is stronger compared to other locations; this also indicates that the problem of stall originates in the vaneless region between the impeller exit and the diffuser inlet. Figure 5-10 Location of the numerical monitor points at the tip clearance gap Figure 5-11 Dynamic pressure variation with time at points in the tip clearance gap 53 Based on the results shown in Fig. 5-11, the vanless region may be the reason for the high pressure fluctuations at the impeller exit; so, a point was selected in the vaneless region to monitor the pressure variation with time, as shown in Fig. 5-12. It may be seen from Fig. 5-12 that the pressure fluctuations develop with time from 250000 Pa to 450000 Pa, and then decrease again to 350000 Pa. This is a part of the surge cycle when the back flow appears at the maximum pressure value in the curve, and this cycle is supposed to be repeated during surge. Figure 5-12 also illustrates velocity vectors for a part of the vaneless region at three moments: before stall initiation, during stall, and during surge. The velocity vectors before stall initiation are uniform without any deviation, but during stall, velocity vectors start to turn from some diffuser passages to the impeller. During surge, the velocity vectors show that there is a complete back flow from the diffuser to the impeller. Pressure (pa) Number of time steps Figure 5-12 Pressure variation with time during stall and surge for a point in the vaneless region 54 During the solution of the operating condition near the surge, the pressure values at the vaneless region between the impeller and diffuser were monitored over time for the pressure fluctuations. Figure 5-13 indicates the pressure values with time (in rotor revolutions) at 12 locations at the vaneless region, and it can be noted that the pressure fluctuates at a small amplitude for some time, and then the amplitude grows gradually over time until it reaches a moment at which pressure starts to fluctuate sharply; this sharp pressure fluctuation indicates that the rotating stall begins at this time and then develops. Pressure at 12 annular positions Time (rotor revolutions) Figure 5-13 Pressure monitoring with time at the vaneless region 5.2 Rotating Stall Simulation for a Compressor with a Vaneless Diffuser and its Comparison to a Compressor with a Vaned Diffuser The numerical model of the compressor with a vaneless diffuser is composed of three main sections: the inlet pipe with a bell mouth attached to its exit, the impeller and the vaneless diffuser with impeller hub cavity, followed by a shaft seal as shown in Fig. 5-14. The straight pipe section, with a length of 1.5 m, was added at the inlet to ensure that the flow is fully developed at the impeller inlet. Also, there is a radial to axial bend added after the vaneless diffuser section in order to redirect the flow in the axial direction. The axial part of the bend was 55 extended to be longer than its actual size to reduce the effect of the outlet boundary condition on the flow distribution inside the vaneless diffuser. Figure 5-14 The numerical model for the compressor with a vaneless diffuser The results compare a compressor with a vaneless diffuser and a compressor with a vaned diffuser, at design condition close to surge, and at surge condition. For the time averaged results section, the sampling time was chosen to be between 0.025 sec and 0.11 sec, which is equivalent to 18 rotor revolutions from the start of the sampling at time, 0.025 sec. The reason for this choice is to ensure the convergence of the solution before sampling. There are several operating conditions used in the simulations in order to test the compressor in a range of conditions, from its design condition to surge. Figure 5-15 shows the location of six different operating conditions: three of the conditions are for the vaneless diffuser case and the other three are for the vaned diffuser case. Operating conditions 1, 2, 3 represent the design condition, close to surge and at surge for the vaneless diffuser case, while the operating conditions 4, 5, 6 represent the design condition, close to surge and at surge for the vaned diffuser case. Full Blade SplitterBlade 56 Figure 5-15 Description of the operating conditions used for the simulations of the vaneless and the vaned diffuser cases 5.2.1 Monitoring the Fluctuating Parameters During Surge (a) Mass flow rate monitoring Mass flow rate monitoring can indicate the degree of flow unsteadiness and reversed flow in the case of stall development. The mass flow rate averaged over the area at the diffuser exit was monitored overtime near surge condition as shown in Fig. 5-16. For the vaneless diffuser case shown in Fig. 5-16 (a), there are strong fluctuations with high amplitudes, and after some time the averaged mass flow rate value changes rapidly, with smaller fluctuations. On the other hand, for the vaned diffuser case shown in Fig. 5-16 (b), similar behaviorappears but with amplitudes lower than those for the vaneless diffuser case; furthermore, each pulse takes more time relative to the vaneless diffuser case. For the vaneless diffuser case, there is no guide for the flow motion at the diffuser as the flow is not forced to move in a specific path; so, the flow moves randomly at the diffuser exit and this random motion results in the rapid mass flow rate fluctuations that occur during stall. For the vaned diffuser case, during stall, there are some diffuser passages affected by stall in the vaneless region, and also, there are some passages not affected by stall. The flow direction at the diffuser exit from all passages is not the same in the circumferential direction and there is a weak flow at the affected passages only; therefore, the 57 mass flow rate fluctuations are relatively small because of the localized reversed flow at the affected diffuser passages during stall. Mass flow rate (lbm/sec) Mass flow rate (lbm/sec) Time (rotor revolutions) Time (rotor revolutions) (a) For the vaneless diffuser case at the operating point 2 (b) For the vaned diffuser case at the operating point 5 Figure 5-16 Area weighted average of the mass flow rate at compressor exit near surge condition for the vaneless and vaned diffuser cases (b) Static pressure monitoring at the area between the impeller and the diffuser The stall initiates inside the impeller due to the formation of the secondary flow and the tip leakage flow, however stall development begins at the diffuser entrance. The chance of stall development occurs at the vaneless region between the impeller and diffuser for the vaned diffuser case, and at any position in the diffuser close to the impeller exit for the vaneless diffuser case. For this reason, it is very important to study the pressure variation of the flow after it leaves the impeller. Static pressure was monitored with time at a location with a radius ratio of 1.1 near the surge condition, as shown in Fig. 5-17, in order to identify the characteristics of the stage before surge. For the vaneless diffuser case, the pressure decreases gradually from 305 kpa to 20 kpa, accompanied with high pressure fluctuations. For the vaned diffuser case, the pressure varies periodically around an average value of 325 kpa for some time, then drops suddenly to 100 kpa, and then begins to rise again gradually. These pressure fluctuations are part of the surge cycle which repeats continuously during surge. Figure 5-18 illustrates the pressure surge cycle for the vaned diffuser case according to the measurements [72]. The numerical results described in Fig. 5-17 (b) are consistent with the measurements shown in Fig. 5-18. The dotted rectangle shown in Fig. 5-18 highlights one of the surge cycles obtained experimentally and it is clear that the pressure curve from numerical results in Fig. 5-17 (b) has nearly the same trend as measurements, except for the moment before 58 the appearance of the pressure drop, where the pressure does not quite reach the 400 kpa level but reaches 380 kpa instead. There is a very interesting difference between the pressure variation during stall for the vaned and vaneless diffuser cases. The main difference is that the pressure decreases gradually with time during stall for the vaneless diffuser case, but decreases suddenly and sharply for the vaned diffuser case, and this is due to the differences in the way the stall cells develop for the two cases. During stall for the case of the vaned diffuser, the flow leaves the impeller with an angle that deviates from the design angle, and then the flow passes through the vaneless region and at some locations, the flow fails to enter the diffuser passages, causing stall cells to be formed at these locations. After some time, some stall cells merge to form a large stall cell, causing a weak flow inside a large number of diffuser passages; at this moment, the reversed flow at the diffuser exit appears as part of a surge cycle with a sharp rise and drop of the pressure value for the case of the vaned diffuser. When stall occurs in the case of the vaneless diffuser, there is a kind of symmetry of pressure values in the tangential direction when the stagnant flow begins to appear; the flow is not guided by vanes as in the vaned diffuser case. This can explain the gradual pressure drop at the beginning of the surge cycle for the case of the vaneless diffuser, as stall develops more uniformly than in the case of the vaned diffuser. Static pressure (Kpa) Static pressure (Kpa) Time (rotor revolutions) Time (rotor revolutions) (a) For the vaneless diffuser case at the operating point 2 (b) For the vaned diffuser case at the operating point 5 Figure 5-17 Static pressure at radius ratio of 1.1 near surge condition for the vaneless and vaned diffuser cases 59 Figure 5-18 Pressure surge cycle for the vaned diffuser case according to measurements [72] (c) Dynamic pressure monitoring at the area between the impeller and the diffuser The dynamic pressure was also monitored with time as shown in Fig. 5-19 and the spectral analysis of the dynamic pressure signal is indicated in Fig. 5-20. It may be noted that the pressure fluctuations in the case of the vaneless diffuser (Fig. 5-19 (a)) are similar to those for the vaned diffuser case (Fig. 5-19 (b)) but the number of pulses in the case of vaneless diffuser is larger in the same time frame. The spectral analysis of the vaneless diffuser case (Fig. 5-20 (a)) shows that there are 2 major peaks and several minor peaks while the spectral analysis of the vaned diffuser case (Fig. 5-20 (b)) indicates that there are more than 2 major peaks and with higher magnitude values. This means that the impact of the dynamic fluctuations in the case of the vaned diffuser is stronger than the impact of the dynamic fluctuations for the vaneless diffuser case, because as the number of frequency pulses increases, this indicates that stall cells are formed and are moving with a specific speed. Dynamic pressure (Kpa) Dynamic pressure (Kpa) Time (rotor revolutions) Time (rotor revolutions) (a) Vaneless diffuser case at the operating point 2 (b) Vaned diffuser case at the operating point 5 Figure 5-19 Dynamic pressure at radius ratio of 1.1 near surge condition for the vaneless and vaned diffuser cases Part of the surge cycle is in a good match with the results shown in Fig. 5-17 (b). 60 Magnite (Pa) Magnite (Pa) Frequency (Hz) Frequency (Hz) (a) Vaneless diffuser case (b) Vaned diffuser case Figure 5-20 Spectral analysis of the convergence history of the dynamic pressure at radius ratio of 1.1 near surge condition for the vaneless and vaned diffuser cases 5.2.2 Velocity Contours Comparison It is also important to study the velocity variation near the impeller entrance as this can indicate the degree of severity of the stalled area formed at the impeller shroud during stall. Figure 5-21 shows the normalized axial velocity contours at a plane located upstream of the impeller inlet for the vaneless and vaned diffuser cases at three different operating conditions: design, near surge, and at surge. The axial velocity is normalized by relating it to the tip velocity at the impeller exit. For the vaneless diffuser case, it may be noted from Fig. 5-21 (a) that the velocity distribution is nearly uniform without stall cells at the design condition. At the near surge condition, the velocity decreases in the radial direction from the hub surface towards the shroud surface, as shown in Fig. 5-21 (b). Also, near surge, there is a very low velocity region formed near the shroud surface and there are some small areas with negative values. Figure 5-21 (c) indicates the situation at surge, in which the low velocity area is increased and the overall axial velocity values are decreased compared to the results of the near surge condition. For the vaned diffuser case, the velocity distribution is uniform at the design condition except for a very small area with a slightly lower velocity value as shown in Fig. 5-21 (d). When the compressor operates near surge, some stall areas are formed randomly with different sizes and intensities and they are localized close to the shroud surface, as illustrated in Fig. 5-21 (e). 61 Figure 5-21 (f) indicates that at surge, the stall areas that appeared earlier are merged, forming larger areas, and these areas extend from the shroud surface toward the hub surface. It may be noted from this comparison that in the case of a compressor with a vaneless diffuser, the stall cells appear near the shroud surface in the form of one continuous layer developing overtime. On the other hand, for the vaned diffuser case, the stall cells appear at random locations and some of them merge together causing the stall area to enlarge and reach the hub surface. For the vaneless diffuser case For the vaned diffuser case va/Utip exit va/Utip exit (a) At design condition (operating point 1) (d) At design condition (operating point 4) (b) Near surge (operating point 2) (e) Near surge (operating point 5) (c) At surge (operating point 3) (f) At surge (operating point 6) Figure 5-21 Axial velocity contours upstream the impeller inlet at different operating conditions for the vaneless and the vaned diffuser cases 62 5.2.3 Pressure Contours Comparison In order to demonstrate the degree of unsteadiness during stall, the Root Mean Square (RMS) of the static pressure at the impeller blades and hub surfaces are shown in Fig. 5-22 for the vaneless and the vaned diffuser cases at the design, near surge and at surge conditions. The higher values of pressure indicate higher unsteadiness or relatively stronger change overtime. For the vaneless diffuser case For the vaned diffuser case bar (a) At design condition (operating point 1) (d) At design condition (operating point 4) (b) Near surge (operating point 2) (e) Near surge (operating point 5) (c) At surge (operating point 3) (f) At surge (operating point 6) Figure 5-22 RMS of the static pressure at the impeller hub and blades surfaces at different operating conditions for the vaneless and the vaned diffuser cases. 63 For the vaneless diffuser case, the high pressure values are concentrated at the impeller blades’ leading edge, and at the blade tips at the design condition as shown in Fig. 5-22 (a). Near surge, the areas of high pressure values are increased, especially at the blade tips near the inlet and exit of the impeller, as indicated in Fig. 5-22 (b). At surge, the area most affected by pressure fluctuations is the impeller exit; the pressure values are slightly higher than those for the near surge condition, as seen in Fig. 5-22 (c). For the vaned diffuser case, at the design operating condition, the pressure values are low at most locations inside the impeller, except near the impeller exit, as shown in Fig. 5-22 (d). Near surge, there are some impeller passages with high pressure values from inlet to exit while there are also passages with relatively low pressure values, as indicated in Fig. 5-22 (e). Figure 5-22 (f) describes the strong pressure instability at the impeller tip exit area at surge condition, with high pressure values formed at the splitter blades’ leading edge up to the impeller exit, where the pressure values are maximized. The comparison of pressure contours reveals that the maximum pressure fluctuations are at the impeller tip exit, and that the pressure distribution is nearly symmetrical for the vaneless diffuser case, while in the case of vaned diffuser, every passage has different pressure distributions. This means that for both the vaned and vaneless cases the most critical location regarding the deformation is at the impeller tip exit because the strong pressure fluctuations cause high stresses on the blades. Consequently, the blade tips may deform inside the tip clearance gap and approach the impeller shroud surface. 5.2.4 Velocity Vectors Comparison Figure 5-23 describes the comparison of the velocity vectors between the vaneless and vaned diffuser cases at a plane located at 98% of the span length, when the compressor is operating at surge condition. For the vaneless diffuser case, the flow inside the diffuser is nearly tangential and turns back toward the impeller at some locations; there is also a very low velocity region near the diffuser exit with stagnant areas, as seen in the left column of Fig. 5-23 in the magnified images. Furthermore, the velocity distribution is identical at impeller exit for all passages. For the vaned diffuser case, the magnified images in Fig. 5-23 in the right column show that some impeller passages are exposed to back flow from the diffuser due to the failure of 64 the flow to enter the diffuser, which is a result of weak flow and low velocity angle at the impeller exit. Also, this back flow causes the stalling of the corresponding diffuser passages and a completely stagnant flow appears inside these passages. Vaneless diffuser case Vaned diffuser case m/s Very low velocity region with back flow Some diffuser passages are stalled Flow direction is nearly tangential Back flow at the diffuser throat Figure 5-23 Velocity vectors at 98% of span for the vaneless and vaned diffuser cases during surge 5.3 Air Injection Simulations The full annulus model was used to simulate the compressor’s behavior during stall. The simulations were performed at an operating condition outside the steady region of the compressor map. Figure 5-24 shows the compressor map (the total pressure ratio versus mass flow rate) with an indication of the selected point for simulations. For the speed curve of 21,789 rpm, which is the design speed for this compressor, there is a point at the end of this speed curve to the left of the surge line corresponding to a mass flow rate of 4.1 kg/s. When the compressor 65 works at a mass flow rate less than this value for the design speed line, the compressor is exposed to stall. Therefore, a point was selected at the left of the surge line corresponding to a mass flow rate of 3.8 kg/s in order to test the compressor behaviorat this point, with and without injection, and at several angles and different injection mass flow rates. Since the numerical model will not converge if we start directly at the point of stall, the model was used first for an operational point close to the surge limit as shown in Fig. 5-24, and the obtained steady solution was then used as an initial value to begin the simulation at the point of stall. Figure 5-24 Location of the point selected for the air injection simulations 5.3.1 Description of the Location and Orientation of the Injected Air In the vaneless region, there is no guide for the motion of the fluid after it leaves the impeller. At low flow rates, there is a probability of reversed flow formation at the diffuser throat due to the random flow motion at the vaneless area. For this reason, the injection of air was selected to be in the vaneless region at the shroud surface in order to minimize the flow disturbances caused by the tip leakage flow. The air was injected at twelve points through holes each with a diameter of 2 mm, as shown in Fig. 5-25 (a). The selection of the number of injection points and the diameter of injection holes for this compressor was based on recommendations clarified in previous studies [22, 46]. The injected air was directed toward the diffuser with an angle Ɵ measured from the tangential direction in the radial-tangential plane as indicated in Fig. 5-25 (b). Several values of Ɵ were selected as follows 0◦, 10◦, 20◦, 30◦, 40◦, 180◦. The injection at Ɵ = 0◦ is known as the forward tangent injection while at Ɵ = 180◦, it is called the reverse tangent 66 injection. The forward tangent injection can increase the stall margin with lower efficiency losses compared to the reverse tangent injection [22]. Therefore, the Ɵ values of 10◦, 20◦, 30◦, 40◦ were chosen in order to find the best injection angle in the range close to the forward injection angle. (a) Injection locations (b) Injection angle Figure 5-25 Description of the injection locations and the injection angle 5.3.2 Clarification of the Air Injection Effect The role of the air injection can be clarified by comparing flow behaviorbetween the case with an injection and the case without an injection. Figures 5-26 (a) and 5-27 (a) show pressure distribution and velocity vectors respectively near the shroud surface during the stall development at a region affected by stall for the case without the use of injections; meanwhile, figures 5-26 (b) and 5-27 (b) show the pressure distribution and velocity vectors respectively for the same region for the case with an injection, at an angle of 20◦ and with a mass flow rate of 1.5%. Figure 5-26 (a) shows that there is a high pressure region at the impeller exit (inside the ellipse shown on the figure); this high pressure region is caused by the back flow coming from the diffuser and impacting the blade as illustrated in Fig. 5-27 (a). The effect of an air injection is clear: the disturbed flow at the vaneless region is corrected as shown in Fig. 5-26 (b) and the velocity vectors’ direction is moved towards the diffuser passage; this leads to normal and identical pressure distribution in the diffuser passages. Direction of rotation θTangential direction Radial direction Injection point67 (a) Without using injection (b) With using injection at angle of 20◦ Figure 5-26 Pressure contours near the shroud surface (a) Without using injection (b) With using injection at angle of 20◦ Figure 5-27 Velocity vectors near the shroud surface 5.3.3 Comparison Between all Injection Cases During stall, the fluid angle at the diffuser inlet is lower than the ideal value, and part of the fluid returns back to the impeller instead of entering the diffuser. The number of diffuser passages with reversed flow was taken as an indication of the degree of success of the air injection to control stall. Figure 5-28 shows the total number of diffuser passages affected by reversed flow during stall for all injection cases at the different injection angles and mass flow 68 rates. It may be concluded from Fig. 5-28 that increasing the injection mass flow rate minimized the total number of reversed flow locations for all injection angles except for Ɵ = 10◦ and Ɵ = 180◦. For the case of Ɵ = 10◦, there is no effect of increasing the injection mass flow rate because in this case, the injected air is directed toward the diffuser’s leading edge point without direct interaction with main flow entering the diffuser. For the case of Ɵ = 180◦, the reversed flow problem worsened as the injection mass flow rate was increased because the injection is administered against the fluid motion; therefore the fluid velocity is excessively decreased. In general, the minimum number of reversed flow passages occurs with an injection angle of Ɵ = 30◦ especially at an injection rate of 1.5% of the design mass flow rate (ṁdesign) where the reversed flow disappeared completely. Figure 5-28 The number of diffuser passages with reversed flow for various angles of injection and injection mass flow rates Figure 5-29 shows the path lines of the injected air for various injection angles and for the minimum and maximum injection mass flow rates at the section of the compressor most affected by stall formation. It may be noted from Fig. 5-29 that all the cases with the minimum injection mass flow rate (0.7% of ṁdesign) have some path lines not directed correctly through the diffuser passage. The worst cases are Ɵ = 0◦, 10◦, 20◦ at an injection rate of 0.7% of ṁdesign,and the case of Ɵ = 180◦ at an injection rate of 1.5% of ṁdesign ,because there are some path lines that transfer tangentially and move randomly towards unexpected diffuser passages or return back to the impeller. 69 Θ = 0◦ at ṁ = 0.7% ṁdesign Θ = 0◦ at ṁ = 1.5% ṁdesign Θ = 10◦ at ṁ = 0.7% ṁdesign Θ = 10◦ at ṁ = 1.5% ṁdesign Θ = 20◦ at ṁ = 0.7% ṁdesign Θ = 20◦ at ṁ = 1.5% ṁdesign Θ = 30◦ at ṁ = 0.7% ṁdesign Θ = 30◦ at ṁ = 1.5% ṁdesign Θ = 40◦ at ṁ = 0.7% ṁdesign Θ = 40◦ at ṁ = 1.5% ṁdesign Θ= 180◦ at ṁ = 0.7% ṁdesign Θ= 180◦ at ṁ = 1.5% ṁdesign Figure 5-29 Flow path lines colored by Mach number value for various angles of injection at the minimum and maximum injection mass flow rate used 70 The deviated path line portions of the worst cases are shown inside the circles in Fig. 5-29. By increasing the injection mass flow rate to 1.5% of ṁdesign instead of 0.7%, the path lines’ directions are corrected, with slight variations, except for the case of an injection angle of 180◦. For the case of injection at the angle of 180◦, the path lines’ directions are deteriorated and shifted in the tangential direction at the injection rate of 1.5% of ṁdesign compared to the injection rate of 0.7% of ṁdesign. This means that increasing the injection mass flow rate is not always beneficial, especially when the injection direction is against the main flow direction. The most successful injection method is at an angle of 30◦ with an injection mass flow rate of 1.5% of ṁdesign because the path lines move through the middle of the diffuser passage without deviation and without the formation of vortices. This means that the injected air in this case is able to modify the fluid velocity and direction entering the diffuser to be close to the design requirements. Based on the comparisons shown in Figures 5-28 and 5-29, the best 6 injection cases among all the 18 injection cases are Ɵ = 0◦, 10◦, 20◦, 30◦, 40◦ at an injection rate of 1.5% of ṁdesign and Ɵ = 180◦ at injection rate of 0.7% of ṁdesign. Also, results indicate that the case of a 30◦ injection angle provides the best results; however, in order to verify this result, the following comparisons were applied to the 6 best injection cases. Figure 5-30 shows one of the paths used for plotting the velocity variation for different injection cases for the specified diffuser passages. This path is composed of two lines; the first line connects the impeller exit with the vaneless region exit passing through the injection point location, and the second line connects the vaneless region exit with the diffuser exit. Figure 5-30 Representation of one of the paths through a diffuser passage used to show results Figure 5-31 indicates the Mach number variation, along the radial distance for different injection mass flow rates and angles of injection (10◦, 20◦, 30◦, 40◦ and 180◦) through the path 71 shown previously in Fig. 5-30. The vaneless region begins at a radius of 0.2175 m and ends at 0.234 m, while the injection point location is located at a radius of 0.227 m, where the velocity peak also occurs due to injection. Injection Rate (0% [no injection])Injection Rate (0.7% of ṁdesign) Injection Rate (1% of ṁdesign) Injection Rate (1.5% of ṁdesign) Figure 5-31 Mach number variation with the radial distance at different injection mass flow rates for injection angles of 10◦, 20◦, 30◦, 40◦ and 180◦ 72 The diffuser passage that was selected to show the effect of increasing injection mass flow rate is different for each case. The diffuser passage selection was based on information in Fig. 5-28, through selecting one of the diffuser passages which stalled at a low injection mass flow rate, but was restored by increasing the injection mass flow rate. The main purpose of the comparison presented in Fig. 5-31 is to show the role of increasing the injection velocity in removing stall for different injection angles. The increase of the injection mass flow rate from 0.7% to 1.5% of ṁdesign resulted in increasing the Mach number at the diffuser entrance (at r = 0.2345 m): from 0.5 to 0.6 in the case of an injection angle of 10◦ (Fig. 5-31 (a)); from 0.55% to 0.63% of ṁdesignin the case of an injection angle of 20◦ (Fig. 5-31 (b)); from 0.61% to 0.75% of ṁdesign in the case of an injection angle of 30◦ (Fig. 5-31 (c)); from 0.5% to 0.7% of ṁdesign in the case of an injection angle of 40◦ (Fig. 5-31 (d)) and from 0.4% to 0.55% of ṁdesign in the case of a reverse tangent injection (Fig. 5-31 (e)). For the case of no injection, the Mach number is 0.25 at the diffuser inlet only, which is smaller than the corresponding value in all injection cases. It may be concluded from Fig. 5-31 (b & c & d) that there is a small velocity peak at the diffuser entrance (at r = 0.2345 m) for an injection rate of 1.5% of ṁdesign; this means that injection at angles of 20◦, 30◦ and 40◦ are the best orientations for modifying the velocity at the diffuser inlet or transferring kinetic energy in the correct direction. The reverse tangent injection case is the worst among all injection cases, because it is based on reducing the tangential velocity in order to correct the velocity angle, but this method is only suitable for working at lower pressure ratios compared with the other injection angles. 5.3.4. Comparison Between the Best Injection Cases In order to identify the best suitable injection method, the best cases for each injection angle were compared. It was found that the best results were obtained at an injection mass flow rate of 1.5% of ṁdesign for the injection angles of 10◦, 20◦, 30◦ and 40◦. It was also shown that for the case of a reverse tangent injection, the best results were obtained at an injection mass flow rate of 0.7% of ṁdesign. 73 (I) Velocity distribution comparison for the best injection cases By combining the best injection cases in one graph (Fig. 5-32), it can be concluded that the worst case is the reverse tangent injection, while the best case is the injection angle of 30◦, with the other three cases rating in-between the two extremes. The differences between these cases can be seen at the end of the vaneless region or at the diffuser inlet, because a jump in velocity magnitude occurs at this location due to air injection. The ordering of the best to the worst cases is as follows: an injection at the angle of 30◦, 20◦,10◦, 40◦, and the reverse tangent injection. Figure 5-32 (a) demonstrates that when the velocity peak at the end of the vaneless region is higher, the entire velocity distribution is shifted to higher values. The fluid velocity angle distribution in the vaneless region for all the selected best injection cases is illustrated in Fig. 5-32 (b) for ease of comparison. The ideal flow occurs when the fluid velocity angle is equal to or close to the diffuser vane angle at the diffuser inlet. Knowing that the diffuser vane angle for this compressor is 13.5◦, it is clear from Fig. 5-32 (b) that the reverse tangent injection is the least effective case because there is a big gap between the fluid and vane angles (about 3.5◦), while the fluid velocity angle for the other remaining injection cases is very close to 13.5◦, the best cases being injections with angles of 30◦ and 40◦. Injection by 1.5 % of ṁdesign & at 10◦Injection by 1.5 % of ṁdesign & at 20◦Injection by 1.5 % of ṁdesign & at 30◦Injection by 1.5 % of ṁdesign & at 40◦Injection by 0.7 % of ṁdesign& at 180◦ (a) Mach number variation (b) Velocity angle variation Figure 5-32 Mach number and velocity angle variations for the best injection cases 74 (II) Pressure contours comparison for the best injection cases At the beginning of stall, low velocity regions are formed at some positions in the vaneless space between the impeller and diffuser near the shroud surface. These low velocity regions are equivalent to high pressure regions at the same time. The high pressure regions are increasing with time during stall development and can be extended to impeller blades close to the vaneless area. This means that the impeller blades are passing through high pressure regions and then through low pressure regions. These pressure fluctuations are strong near the impeller exit and can break the impeller blades. Figure 5-33 shows the pressure distribution at a plane close to the impeller shroud surface for the best five cases of injection selected earlier, along with the case of no injection after 5 revolutions from the start of instabilities and pressure fluctuations. The following key points are concluded from Fig. 5-33: - At an injection angle of 10◦, two slightly high pressure regions occur at the impeller exit. There are no high pressure regions occurring at an injection angle of 20◦ and 30◦. The case with an injection angle of 30◦ is the best in achieving the compressor stability because the pressure distribution is more uniform radially and the pressure values are lower than those at an injection angle of 20◦. - By increasing the angle of injection to 40◦, there is only one high pressure region, and for the case of a reverse tangent injection, there are two large regions of high pressure. - The worst case scenario is the case of no injection, where the impeller blades are exposed to extreme and deep pressure fluctuations compared to cases in which injections are used. 75 For injection angle = 10 ◦ For injection angle = 20 ◦ For injection angle = 30 ◦ For injection angle = 40 ◦ For reverse tangent injection For the case of no injection Figure 5-33 Pressure contours (in Pascals) at impeller shroud for the best injection cases (III) Forces acting on the impeller for the best injection cases The net forces acting on the compressor rotator vanes were calculated by the integration of pressure on the suction and pressure sides, and the tip surface of the impeller blades. The forces are calculated in the x, y and axial directions. These directions are shown in Fig. 5-34. 76 Figure 5-34 The CFD numerical model of the NASA CC3 compressor Figures 5-35 (a), 5-35 (b) and 5-35 (c) show the variation of forces overtime (in rotor revolutions) in the axial, x and y directions respectively for the best injection cases at different injection angles, and also for the case with no injection. The reference time (t=0) is the time at the beginning of the appearance of high pressure fluctuations in the vaneless region in the case of no injection. Figure 5-35 also shows that for the case of no injection, the force magnitudes are large compared to those for the other cases of injection, while the smallest force values occur with the injection at angles of 20◦ and 30◦. Figure 5-35 (a) indicates that the injections at angles of 20◦ and 30◦ are the best cases for the stabilization of the compressor because the force distribution is more uniform with small fluctuations. The case of a reverse tangent injection, an injection at an angle of 10◦, and the case of no injection are the worst cases because that the force increases gradually with time, especially for the case with no injection, where the largest force values occur with strong fluctuations. Figure 5-35 (b) indicates that the best case is the injection with an angle of 30◦ because the force distribution has the lowest fluctuation amplitude, nearly fluctuating around zero. Also, the case with an injection angle of 20◦ shows that the force fluctuates around some value which is not equal to zero. For the case of no injection, the force changes randomly and reaches extreme values as shown in Fig.5-35 (c). All the other cases show that the force increases with time except for the Bend Full impeller blade Splitter impeller blade Diffuserx-directiony-directionAxial direction 77 case of an injection angle of 30◦ at which there are fluctuations between 0 &100 N, and the average force is not increasing with time as in the other cases. In summary, it may be concluded that the best stabilized case is the injection with an angle of 30◦ because the force distribution is nearly uniform, with the smallest fluctuations compared to the other cases in which the average force increases gradually over time. (a) In the axial direction (b) In x- direction (c) In y- directionFigure 5-35 Variation of the force acting on blades with time for the best injection cases and for the case without using injection 78 5.3.5. Effect of Changing the Injection Angle in the Radial-Axial Plane The variation of the injection angle Ɵ at the radial-tangential plane was studied in detail as discussed above. Another factor was taken into consideration, which is studying the effect of changing the injection angle at the radial-axial plane. Figure 5-36 indicates an explanation of the injection angle α defined at the radial-axial plane where the fluid is injected at a point at the shroud surface and is directed toward the diffuser inlet with an angle α measured from the radial direction. All previous results were calculated at α = 5◦ and the best obtained results were achieved for the case of Ɵ = 30◦. Figure 5-36 The angle of injection α at the radial-axial plane The following results are for the case of Ɵ = 30◦ at two other values for α (10◦ and 15◦) in comparison with the previous results obtained at α = 5◦. Figure 5-37 (a, b and c) shows the variation of pressure with time for the injection cases of α = 5◦, 10◦, 15◦ at Ɵ = 30◦ for points located at the vaneless region at 90% of the span length. It may be seen from Fig. 5-37 (a) that the pressure fluctuations for all points are nearly identical, and also the minimum and maximum pressure values are nearly constant. For the case of α = 10◦ shown in Fig. 5-37 (b), the pressure fluctuations are close to those of the case of α = 5◦, but there are two signals with a sudden change of pressure values. The amplitudes of fluctuation of these two signals are higher than those for the other signals; this means that at the location of these two points in the vaneless region, a stall area is initiated. For the case of α = 15◦ as shown in Fig. 5-37 (c), the pressure fluctuations for all signals are random and with relatively high amplitude; this ensures that the 79 injection at α = 15◦ fails to keep the fluid stabilized in the vaneless region compared to other cases. Based on the obtained results, the injection of air with Ɵ = 30◦ and α = 5◦ at an injection mass flow rate of 1.5% of ṁdesign achieved compressor stability at an operating condition in the unstable zone corresponding to a mass flow rate of 3.8 kg/s. (a) For α = 5◦ (b) For α = 10◦ (c) For α = 15◦ Figure 5-37 Pressure variation at points in the vaneless region with time at Ɵ = 30◦ for α = 10◦, 15◦, 20◦ 80 5.4 Casing Treatment Simulations This section is divided into two segments; the first segment (5.4.1) discusses the simulations of various casing treatment configurations whereas the second segment (5.4.2) is about a specific type of casing treatment: casing grooves, also known as casing slots. The aim of the first segment is to compare the efficiency of different casing treatment methods for reducing stall. The main goal of the second segment is to study the effects of changing some of the parameters on the performance of the casing grooves. 5.4.1 Simulation of Different Casing Treatment Configurations There are four numerical models used in the simulations of this segment. The first model represents the compressor without modifications in its casing. The second model contains two types of casing treatments, while the third model is identical to the second model but enhanced with an external air injection. The fourth model contains one type of casing treatment. The first model is shown in Fig. 5-38 (a); there is no casing treatment in this model. It is composed of four main components: the inlet guide duct, the impeller, the diffuser and the shroud thin layer. The thin layer surface is smooth for the first model, but for the other three models, this layer is modified according to the casing treatment type. The function of the shroud thin layer is to capture the tip leakage flow and its thickness is 0.1 mm, which is half of the tip clearance value. The other half of the tip clearance belongs to the impeller rotating zone which slides on the thin layer fixed zone. The second model combines the holed casing treatment and the casing slots (Fig. 5-38 (b)); this model is called case 1 in this study. There are 15 tubes connecting the impeller casing with the inlet guide duct. Each tube has a diameter of 3 mm and this diameter was selected based on a recommendation from previous studies [20, 44]. There is a slot in the casing and the dimensions of this slot are assumed to be 4 mm in width and 2 mm in height. The third model is the same as the second model, but the slot is modified to have 15 injection holes through its upper wall. Figure 5-38 (c) shows the configuration of the slot with the air injection positions. This model is called case 2 in this study. The air is injected through these holes in the tangential direction at an injection mass flow rate of 1% of the inlet design mass flow rate, and this value was assumed because it is the middle of the range used in several previous studies [3, 22, 70, 78, 79]. The fourth model contains three casing slots with the same dimensions as described in the 81 first and the second models; the fourth model is called case 3 in this study. Figure 5-38 (d) describes the shape of the fourth model. The position of slot 1 is the same as in previous models, the position of slot 2 is at the leading edge of the splitter blades, and the position of slot 3 was selected so that that slot 2 is located equidistant between slots 1 and 3. (a) The numerical model without casing treatment (case 0) (a) The combined model of the holed casing and the casing slot (case 1) (c) The combined model of the holed casing and the casing slot with air injection (case 2)(d) The numerical model with casing slots only (case 3) Figure 5-38 Numerical models for three different casing treatment cases and for the case without casing treatment The problem of stall starts near the impeller shroud surface due to the formation of the tip leakage flow. So, it is important to study the velocity contours at the blades tip in order to determine the stall areas and its propagation level. Figure 5-39 (a to d) illustrates the velocity 82 contours averaged over a 10 rotor revolution time period starting from the occurrence of stall at 0.99 of the span length for the case with no casing treatment, case 1, case 2 and case 3, respectively. It can be noted from Fig. 5-39 (a) that there are large areas with relatively low velocity values or stagnation regions which cause a blockage of the flow inside the impeller passages. Figure 5-39 (b) shows that the velocity values are high at the slot position due to the circulation of the tip leakage flow inside the slot, and that the low velocity regions are concentrated at the impeller blades’ leading edge. The overall velocity distribution is enhanced compared to the case with no casing treatment, especially at the area around the slot position and also at the area between the slot position and the impeller exit. m/s (a) For case 0 (b) For case 1 (c) For case 2 (d) For case 3 Figure 5-39 Velocity contours at 0.99 of span for the four different cases 83 The area marked by a circle in Fig. 5-39 (b) is explained in a separate figure above Fig. 5-39 (b); it shows the effect of the holed casing treatment and the velocity vectors inside one of the tubes at the casing. The direction of the flow inside the tube is from the impeller casing towards the impeller inlet, which means that the tip leakage flow is sucked into the tube due to the relatively high pressure area formed at the shroud surface. Figure 5-39 (c) indicates similar conclusions made for Fig. 5-39 (b) but with higher velocity values at the slot position, and this is due to the use of air injections at some points in the slot wall. It can be noted that the velocity distribution is more uniform compared to Fig. 5-39 (b) especially near the impeller trailing edge, and also the velocity values are slightly higher before and around the slot position. Figure 5-39 (d) shows the influence of using several slots distributed along the casing length as the low velocity regions disappear around each slot position. Also, the velocity is higher at the leading edge compared to all other cases discussed before. Casing slots increase kinetic energy due to the circulation effect of the flow at the areas of stagnation. Figures 5-40, 5-41 indicate the kinetic energy gain due to using casing slots by comparing the casing treatment cases with the case without casing treatment. Figure 5-40 shows the average Mach number at the position of the first slot and its variation with time for all cases. It can be noted from Fig. 5-40 that the Mach number for cases 1, 2 and 3 is higher than that for case 4, without casing treatment. Also, the Mach number distribution for cases 1 and 3 are close to each other but case 2 achieved a significant increase in the kinetic energy because of the use of an air injection. Figures 5-41 (a), 5-41 (b) show the variation of the average Mach number over time at the position of the second slot and the third slot, respectively. Figure 5-41 (a) proved again that the kinetic energy is increased at the slot position but this increase is less than that for the first slot. Also, for the third slot (Fig. 5-41 (b)), there is nearly no benefit from the casing slot as the two curves shown are nearly identical. Therefore, it can be concluded that as the position of the slot approaches the impeller inlet, the kinetic energy gained from it decreases. 84 The isentropic efficiency also varies with time during stall. Figure 5-42 indicates the variation of the total to total isentropic efficiency with time during stall. It can be noted from Fig. 5-42 that all casing treatment cases have an efficiency value lower than the case without casing treatment and this is due to the interaction between the main flow in the impeller passages and the secondary flow in the casing slots. Also, it can be seen that case 3 is the best among all of the casing treatment cases. Figure 5-40 Average Mach number variation with time at the position of the first slot (a) At the position of the second slot (b) At the position of the third slot Figure 5-41 Average Mach number variation with time 0.350.370.390.410.430.450.471 2 3 4 5 6 7 8 9 10Average Mach number at the position of the first slotTime (revolutions)For case 1For case 2For case 3For the case of no casing treatment00.050.10.150.20.250.31 2 3 4 5 6 7 8 9 10Average Mach number at the position of the second slotTime (revolutions)For case 3For the case of no casing treatment00.050.10.150.20.251 2 3 4 5 6 7 8 9 10Average Mach number at the position of the third slotTime (revolutions)For case 3For the case of no casing treatment85 Figure 5-42 Isentropic total to total efficiency variation with time during stall for all cases 5.4.2 Simulation of the Casing Grooves with the Variation of their Parameters This section is composed of four parts: the first part discusses the simulations of groove aspect ratio variation; the second part presents the simulations of the effect of changing the groove location; the third part presents the simulations for changing the number of grooves; and, the last part introduces a comparison between the groove flow performance at the best efficiency point and that for the stall operating condition. The stall operating condition is located at the surge line corresponding to mass flow rate of 4.1 kg/s. The numerical model used for the groove aspect ratio simulations is shown in Fig. 5-43. It is indicated in Fig. 5-43 that there is only one groove at the leading edge of the impeller blades and the cross section of this groove has a width (W) and height (H). The width of the groove was assumed to be constant in all cases and equals 5 mm, while the groove height was varied, with assumed values of 1, 3, 5, 7 and 9 mm in order to generate grooves with 5 different aspect ratio (H/W) values equal to 0.2, 0.6, 1, 1.4 and 1.8. For the simulations regarding the effect of changing groove location, there are three numerical models used. The first and second models are for the groove locations 1 & 2 where the groove is located at the impeller full blades leading edge and at the impeller splitter blades leading edge, as seen in Fig. 5-44 (a) and Fig. 5-44 (b), respectively. For the third model (Fig. 5-0.670.680.690.70.710.720.731 2 3 4 5 6 7 8 9 10Total to total efficiencyTime (revolutions)For case 1For case 2For case 3For the case of no casing treatment86 44 (c)), the groove is located at location 3 after the splitter blade leading edge, and it is located in a way that groove location 2 is in the middle of the meridional distance between groove locations 1 & 3. The numerical models used for the simulations of the effect of changing the number of grooves, are shown in Fig. 5-45. Figures 5-45 (a) to (d) show the numerical model with 1 groove, 3 grooves, 5 grooves, and 7 grooves, respectively. It is apparent in Fig. 5-45 that the grooves locations are concentrated near the impeller inlet more than at the impeller exit; this was chosen based on the results extracted concerning the best groove location, and will be discussed in detail. WxH ≡ 5 x 1 mm WxH ≡ 5 x 3 mm WxH ≡ 5 x 5 mm WxH ≡ 5 x 7 mm WxH ≡ 5 x 9 mm Figure 5-43 The numerical model used for the simulations of the groove aspect ratio variation H W Groove 87 (a) location 1 (b) location 2 (c) location 3 Figure 5-44 The numerical model used for the simulations of the groove location variation (a) with 1 groove (b) with 3 grooves (c) with 5 grooves (d) with 7 grooves Figure 5-45 The numerical models used for the simulations of changing the number of grooves 88 (a) Results of the simulations for changing groove cross section aspect ratio Figures 5-46 (a) to (e) illustrate the velocity contours at a plane perpendicular to the axial direction and passing through the groove, with clarification of the velocity vectors inside the groove area for the groove cross section aspect ratios of 0.2, 0.6, 1, 1.4 and 1.8. Fig. 5-46 contains other explanatory figures marked by arrows in order to show the velocity vectors at the grooves cross section at specified locations; these locations are at the start of the reinjected flow coming from the grooves through the impeller passages. In general, the blade pushes the flow inside the groove and this flow is re-injected again from the groove after the blade passes. It may be noted from Fig. 5-46 (a) that the reinjected flow from the groove is weak because the height of the groove is relatively too short to allow the flow to circulate inside the groove properly as shown in right column in Fig. 5-46 (a). For the case of a groove aspect ratio of 0.6 (Fig. 5-46 (b)), the flow has more space inside the groove to move and circulate; as a result, the reinjected flow is more uniform when entering the impeller passage compared to the case in Fig. 5-46 (a); however, the reinjected flow velocity is still small. For the case of a groove aspect ratio of 1 (Fig. 5-46 (c)), the reinjected flow velocity becomes higher, as is apparent when comparing the images on the right side of Fig. 5-46 (c) and Fig. 5-46 (b). When the groove aspect ratio is increased to be more than 1, the reinjected flow velocity increases slightly, and this may be concluded by comparing the images to the right in Fig. 5-46 (d) and Fig. 5-46 (e) with the other previous aspect ratio cases. However, at the same time, the reinjected flow effective area increases in the circumferential direction as is visible when comparing the images on the left of Fig. 5-46 (d) and Fig. 5-46 (e) with previous aspect ratio cases. 89 (a) For the groove aspect ratio of 0.2 (b) For the groove aspect ratio of 0.6 (c) For the groove aspect ratio of 1 (d) For the groove aspect ratio of 1.4 (e) For the groove aspect ratio of 1.8 Figure 5-46 Velocity contours at a plane perpendicular to the axial direction with the clarification of the velocity vectors inside the groove for various groove aspect ratio values Direction of rotation 90 Figure 5-47 shows the plotting of the radial velocity on the vertical axis against the circumferential angle on the horizontal axis for a line that passes circumferentially at the impeller shroud level through the groove for 3 cases with groove aspect ratios of 0.2, 1 and 1.8, which represent the minimum, average and maximum studied values, respectively. The circumferential angle shown in Fig. 5-47 is for a quarter of the circumferential length because it was found to be enough to make the comparison as the complimentary part of the circumferential length has similar trends and values. Figure 5-47 Radial velocity variation with the circumferential angle for 3 different groove aspect ratio values The positive radial velocity values represent the flow entering the groove coming from the area at the blade pressure side. The negative radial velocity values represent the reinjected flow, which leaves the groove at the area near the blade suction side. It may be noted from Fig. 5-47 that the negative velocity values for the case with a groove aspect ratio of 0.2 is relatively low compared to the other cases. For the cases with groove aspect ratios of 1 and 1.8, the maximum negative velocity values are very close to each other. But the area under the curve for the negative velocity values in the case of a groove aspect ratio of 1.8 is higher than that for the case with a groove aspect ratio of 1. In conclusion, Fig. 5-47 demonstrates that when the groove aspect ratio is increased to be more than 1, the reinjected flow velocity remains nearly the same, but the reinjected flow effective area is increased. But, when the groove aspect ratio is decreased 91 to be less than 1, both the reinjected flow velocity and the reinjected flow effective area are decreased. (b) Results of the simulations of changing groove location Figures 5-48 (a), (b) and (c) indicate the Mach number values at a plane located at 99% of the span length of the impeller for the cases with groove locations 1, 2 and 3, respectively. It may be noted from Fig. 5-48 that for the case with groove location 1, the low Mach number values at the entrance are limited to small areas around the blade surface. For the case with groove location 2, the low Mach number area is expanded in the circumferential direction at the impeller leading edge but is minimized in the axial direction compared to the case of groove location 1. For the case with groove location 3, the low Mach number area covers a wide region at the impeller entrance and connects the leading points of the blades, but the area around the groove has high Mach number values compared to the other 2 cases in Fig. 5-48 (a) and Fig. 5-48 (b). Based on the results shown in Fig. 5-48, it may be said that groove location 1 has, relatively, the best results; this is because the flow separation at the inlet is controlled. Meanwhile, for groove location 3, the flow is improved only around the groove location but there is no control of the flow at the impeller entrance. (a) For the groove location 1 (b) For the groove location 2 (c) For the groove location 3 Figure 5-48 Mach number values at 99% of the impeller span for 3 different groove locations Groove 1 Groove 2 Groove 392 In order to clarify the comparison of groove locations, the static pressure value is plotted with the circumferential angle at a line that extends circumferentially below the shroud surface at 99% of the span length at the impeller blades’ leading edge, for the three groove location cases, as shown in Fig. 5-49. Also, Fig. 5-49 demonstrates that all curves have pressure fluctuations and that this is normal because there is a pressure variation between the blade pressure side and the blade suction side. The lowest pressure values are for the case with groove location 1; groove location 1 provides the lowest blade loading at the impeller inlet area, compared to the other 2 cases, as shown in Fig. 5-49. Figure 5-49 Static pressure variation in the circumferential direction at the blades leading edge at 99% of span (c) Results of the simulations in which the number of grooves are varied When the number of grooves is increased, velocity distribution and direction may differ. In order to illustrate this distribution, the stream lines colored by velocity magnitude at 98% of the span length of the impeller are shown in Fig. 5-50((a) to (e)) for the cases including no grooves, 1 groove, 3 grooves, 5 grooves, and 7 grooves. Figure 5-50 (a) indicates that there are 93 strong flow disturbances, especially at the impeller blades’ leading edge, and that the stream lines deviate from the blades in many locations. If the results of the case without grooves were taken as a reference for the other results, the effect of increasing the number of grooves may be realized clearly. Figure 5-50 (b) shows that the stream lines distribution for the case of 1 groove is not improved compared to that for the case without grooves in Fig. 5-50 (a). By using 3 grooves (Fig. 5-50 (c)), the velocity values are increased and the stream lines distribution at the inlet is improved, but flow separations still occur at the middle of the meridional path at the splitter leading edge. The blades’ flow separations are shifted further towards the impeller exit in the case of using 5 grooves as indicated in Fig. 5-50 (d). Figure 5-50 (e) shows that the stall problem is nearly resolved and the flow appears to be uniform with minor disturbances when using 7 grooves. The results shown in Fig. 5-50 give the indication that the location of the flow separations can be shifted towards the impeller exit by increasing the number of grooves. Therefore, the importance of groove location must be taken into consideration because results prove that the stall can be reduced when the groove is located near the blades’ leading edge. (a) For the case of smooth casing (b) For the case of 1 groove (c) For the case of 3 grooves (d) For the case of 5 grooves (e) For the case of 7 grooves Figure 5-50 Velocity stream lines at 98% of span for cases with different numbers of grooves 94 It is important to evaluate the flow behaviorat the area between the hub and shroud surfaces in order to understand the depth of the stall areas. The average velocity contours at the meridional plane of the impeller are shown in Fig. 5-51((a) to (e)) for the cases with no grooves, 1 groove, 3 grooves, 5 grooves and 7 grooves. The worst results are for the case without grooves (Fig. 5-51 (a)) where the stall areas with low velocity values appear near the hub and shroud surfaces. By using 1 groove (Fig. 5-51 (b)), there is still a large stall area at the tip, but the overall velocity values are increased compared to the case without grooves. By tracking the results inside the oval shape added to each image in Fig. 5-51, it can be seen that by increasing the number of grooves, the stall area at the tip region is decreased in size and is shifted towards the impeller blades’ trailing edge at the same time. (a) For the case of smooth casing (b) For the case of 1 groove (c) For the case of 3 grooves (d) For the case of 5 grooves (e) For the case of 7 grooves Figure 5-51 Average velocity contours at the impeller meridional plane for cases with different numbers of grooves Based on the results from Fig. 5-50 and Fig. 5-51, the cases with 5 and 7 grooves have the best stream lines and velocity distributions compared to the other cases. Figure 5-52 presents the average pressure distribution along the meridional direction in the impeller for the case of smooth casing, along with cases in which 5 and 7 grooves are implemented. It may be noted from Fig. 5-52 that the highest pressure values are for the case of the smooth casing, without the use of any grooves. Also, it may be seen from Fig. 5-52 that the pressure values for the case of 7 grooves are slightly higher than those for the case of 5 grooves at up to 30% of the meridional distance, but the contrary occurs in the range of 30% to 100% of the meridional distance. 95 Figure 5-52 Pressure variation with the normalized meridional distance for the cases of smooth casing and for the cases with 5 and 7 grooves Changing the number of grooves does not only affect the velocity and pressure distributions, but also affects the overall isentropic efficiency of the compressor. Figure 5-53 describes the total to total isentropic efficiency at an operating point at the surge limit corresponding to a mass flow rate of 4.1 kg/s for the case of smooth casing and also for all other cases in which grooves are used. There is nearly no difference in efficiency values between the case of smooth casing and the case of using 1 groove. This indicates that the effect of casing grooves on efficiency losses is nearly negligible when using 1 groove only, as shown in Fig. 5-53. It may be noted from Fig. 5-53 that the efficiency starts to decrease gradually when using more than 1 groove. Therefore, it is clear that increasing the number of casing grooves improves the flow distribution near the shroud surface and minimizes the tip leakage flow, but this is accompanied with some losses in the isentropic efficiency. 96 Figure 5-53 Total to total isentropic efficiency for cases with different numbers of grooves (d) Groove performance comparison at different operating conditions The groove performance depends on the compressor operating condition because the groove flow is generated based on the magnitude and direction of the main flow velocity in the impeller passages. Figures 5-54 & 5-55 describe the velocity vectors near the shroud surface and also inside the grooves at the best efficiency operating condition, and at the stall operating condition, respectively. Figure 5-54 Velocity vectors inside the grooves and at a plane near the shroud surface at the best efficiency operating condition Groove flow Main flow Groove flow interferes with the main flow with weak interaction.97 Figure 5-55 Velocity vectors inside the grooves and at a plane near the shroud surface at the stall operating condition At the best efficiency condition (Fig. 5-54), the flow in the groove is relatively weaker than the main flow, and the interaction between the main flow and the groove flow does not affect the main flow direction or magnitude. This interaction depends on the location of the groove; and a slight interaction may occur if the groove is placed at the inlet. Figure 5-55 shows that the flow distribution is different from that at the best efficiency condition and there is a secondary flow appearing: the tip leakage flow. The tip leakage flow is a part of the main flow, and moves from one passage to another through the tip clearance, against the groove flow direction. The magnified views of Fig. 5-55 show two different scenarios that result from the interference between the tip leakage flow and the groove wall. The magnified Main flow Tip leakage flow Groove flowGroove flow Tip leakage flow The flow regenerated correctly inside the impeller passage after strong interaction. Vortex is formed because the tip leakage flow is slightly stronger than the groove flow. 98 view at the right of Fig. 5-55 shows that the tip leakage flow velocity is slightly higher than the groove wall velocity, resulting in the formation of some vortices with tangential flow. In this situation, the reversed flow is prevented, but also, the resulting flow is not in the correct direction. The magnified image at the bottom of Fig. 5-55 shows that the groove flow is stronger than the tip leakage flow, and after the interaction, the groove flow was able to overcome the tip leakage flow, resulting in a flow that is in the correct direction but weaker than the ideal flow. So, the compressor flow stability depends on how much flow escapes from the main passages and forms the tip leakage flow, and also on the groove reinjected flow, which varies according to variations in groove geometry and location. 5.5 Structural Analysis of the NASA CC3 Compressor The structural analysis section presents the results extracted from the mechanical model, including deflections and stresses which were obtained after transferring the aerodynamic load from the CFD model. The results in this section are presented in the form of comparisons between the original designed compressor (the uncontrolled compressor) and the compressor using air injection (the controlled compressor). 5.5.1 Stress Variation Near Surge for the Controlled and the Uncontrolled Compressors The maximum principal stress and the equivalent stress distributions at the blade surface at the steady state condition were discussed in the validation results (section 4.3). During stall, these distributions vary strongly over time. For this reason, the peak values of these distributions were traced during the transition from stall to surge. Figures 5-56 (a) and 5-56 (b) show the variation of the peak value of maximum principal stress with time (in rotor revolutions) for the full blade and the splitter blade, respectively. Figures 5-57 (a) and 5-57 (b) indicate the variation of the peak value of the equivalent stress with time (in rotor revolutions) for the full blade and the splitter blade, respectively. The time frame for Fig. 5-56 and Fig. 5-57 begins at a moment during stall (the reference time) and spans up to the time corresponding to the peak of the first surge cycle. It may be seen from Fig. 5-56 (a) that 99 the maximum principal stress for the full blade, for the two cases, fluctuates around two convergent mean values; however, there is some difference near the end of the curve as the final value of the uncontrolled case reaches 370 Mpa, compared with the final value of 250 Mpa for the controlled case. Fig. 5-56 (b) shows that the stress distribution of the splitter blade for the uncontrolled case reaches relatively high values up to 500 Mpa, but the controlled case maintains the same range of stress values for the full blade. This means that the controlled case is stable because both of the full and splitter blades have nearly the same stress range of variations. (a) For the full blade pressure side (b) For the splitter blade pressure side Figure 5-56 The variation of the peak value of the maximum principal stress with time for the impeller blade The maximum equivalent stress variation for the full blade in the uncontrolled case is slightly higher than for the controlled case (Fig. 5-57 (a)). However, for the splitter blade, there is a sudden variation at surge after 6 rotor revolutions for the uncontrolled case, as shown in Fig. 5-57 (b). Finally, it may be noted that for the uncontrolled case, there are stress fluctuations during stall and then the stress increases sharply when approaching surge. For the controlled case, the stress changes slightly with time and is more stable. 100 (a) For the full blade pressure side (b) For the splitter blade pressure side Figure 5-57 The variation of the peak value of equivalent stress with time for the impeller blade 5.5.2 Total Force and Moment Variations Near Surge for the Controlled and the Uncontrolled Compressors The total force and moment were computed for the uncontrolled and controlled compressor cases in x, y and z directions. The z-direction represents the axial flow direction while the x-y plane is perpendicular to the z-direction. Figure 5-58 ((a) to (e)) presents the variation of Fx, Mx, Fy, My, Fz and Mz respectively with time for both of the controlled and the uncontrolled compressor cases. The time frame shown in Fig. 5-58 is the same as explained for Fig. 5-57, where F stands for the force and M represents the moment. There are strong fluctuations for Fx with time for the case of the uncontrolled compressor (Fig. 5-58 (a)) and also for Fy (Fig. 5-58 (c)); this is due to the non-uniform pressure distribution during stall as a result of the formation of stall cells at some impeller passages. Figure 5-58 (e) shows that Fz for the uncontrolled compressor case increases gradually with time due to the expansion of the stagnation areas at the impeller exit at the curved part of the hub surface. On the other hand, the force fluctuates around some mean values which are nearly constant with time for the controlled compressor case. Also, the mean fluctuation value is nearly zero for the Fx variation, which is a good indicator that the compressor is balanced. The moment variation is a reflection of the fluctuations of the force magnitude and direction. For this reason, the Mx and My variations are close to zero for the controlled compressor case but they have relatively high values for the 101 uncontrolled compressor case. Figure 5-58 (f) shows that there is no significant difference between the two cases for the variation of Mz. In general, it may be concluded that the force and moment values are lower for the case with the use of air injection, compared to the uncontrolled compressor results. (a) Fx variation with time (b) Mx variation with time (c) Fy variation with time (d) My variation with time (e) Fz variation with time (f) Mz variation with time Figure 5-58 The variation of forces and moments with time for both of the controlled and the uncontrolled cases 102 5.5.3 Deflection of the Blades During Stall for the Controlled and the Uncontrolled Compressors The most important results concern the blades’ deflection especially at the blades’ tip zone, because the compressor can be damaged if the deflection value is larger than the tip clearance gap distance. Figure 5-59 (a) shows the full blades’ tip profile with maximum deformation which occurs at the peak pressure of the surge cycle for the uncontrolled compressor case. The base surface profile (before deformation) and the shroud surface profile are also shown in Fig. 5-59 (a). The tip clearance value is relatively very small compared to the complete axial distance, Z, and the radial distance, R. For this reason, Fig. 5-59 (a) does not clearly illustrate the three profiles. Three zones were selected to illustrate the deformation more clearly: zone 1 (Fig. 5-59 (b)) is at the impeller exit, while zone 2 (Fig. 5-59 (c)), and zone 3 (Fig. 5-59 (d)) are at the transition from the axial to the radial profile. In fact, the entire profile was checked for the deformation values, but these three zones were found to be the most critical areas where the deformed blade surface approaches the shroud surface. At zone 1, the deformed full blade surface is close to the shroud surface but it does not touch it, and as seen in Fig. 5-59 (b), there is a deformation in the radial direction; this deformation is not dangerous because it actually deforms freely in the vaneless region, radially. At zones 2 and 3, the deformed full blade surface approaches the shroud surface at a dangerously close proximity, and touches the shroud surface at some points, especially in zone 2 at the z-range from 0.1015 m to 0.103 m, as seen in Fig. 5-59 (c). It was also found that for the splitter blade, the z-range from 0.1017 m to 0.1028 m is a dangerous deflection zone. The reason zone 2 has been identified as the most critical deflection area is that in this region, the radial and axial deflection components are combined together and the resultant component direction is towards the shroud surface. For the inlet impeller region, the profile is mostly axial, so the axial deformation is not effective. For the impeller exit region, the profile is mostly radial, and only the axial deformation is effective. 103 Figure 5-59 The maximum full tip blade surface deformation profile for the case of an uncontrolled compressor For the case of a controlled compressor, the results are much better than those for the uncontrolled compressor case in terms of the maximum deformation values. Figure 5-60 presents the maximum full tip blade surface deformation profile for the case of a controlled compressor at zone 2. By comparing the blade profile in Fig. 5-60 to that in Fig. 5-59 (c), it may be noted that for the controlled compressor case, there is no contact between the blade surface and the shroud 104 surface. This means that using the air injection control was beneficial in decreasing the pressure fluctuations, and hence decreasing blade deformation during stall. Figure 5-60 The maximum full tip blade surface deformation profile for the case of controlled compressor using air injection at zone 2 5.5.4 Effect of Variation in the Bearing Stiffness on the Deformation of the Blades Tip Surfaces During Stall There is a relationship between the blades’ deformation and the compressor shaft fixation method. This section discusses how the deformation of blades is affected by the shaft stability governed by bearings. Figure 5-61 (a) describes the two locations of the bearing while Fig. 5-61 (b) indicates the way that each bearing is numerically represented. Each bearing was assumed to be represented by 12 springs connecting the shaft to reference points which have zero displacement. There are six different values for the spring stiffness assumed in order to cover the range of relatively low to high stiffness values. The ball bearing radial stiffness varies from 8x107 N/m to 108 N/m for the radial load range between 1000 N and 1500 N according to measurements [80] and numerical results [81]. The six values for the total stiffness were selected especially to cover this range, and at the same time, to extend slightly out of this range in order to study the effect of using very low and very high stiffness values. The range of stiffness values 105 selected are 4x107 N/m, 6x107 N/m, 8x107 N/m, 10x107 N/m, 12x107 N/m and 14x107 N/m. The stiffness of each spring was assumed to be equal to the value of the total stiffness divided by the number of springs. The values of total stiffness and spring stiffness are shown in table 5-1. (a) location of the bearings (b) representation of each bearing by springs Figure 5-61 The location of the compressor shaft bearings and their numerical representation Table 5-1 Stiffness values Total stiffness (N/m) 4x107 6x107 8x107 10x107 12x107 14x107 Spring stiffness (N/m) 1.666x106 [ks1] 2.5x106 [ks2]3.333x106 [ks3]4.166x106 [ks4]5x106 [ks5] 5.833x106 [ks6] Figure 5-62 shows the deformed blades tip surface profile at the peak of surge at zone 2 for various bearing stiffness values. It may be noted from Fig. 5-62 (a) that the blade profile touches the shroud surface profile for nearly half of the z-range in zone 2 when using the k1 stiffness value. There are also many points of contact for the other zones, and this indicates that using a ks1 value is not safe for the blades unless the current tip clearance is increased to avoid blade damage. By increasing the stiffness to the value of ks2, the blade deformation profile is slightly shifted but still very close around the position of z = 0.102 m as shown in Fig. 5-62 (b). At the k3 value (Fig. 5-62 (c)), the blades profile comes close to the mid position between the base and shroud surface profiles. Based on this result, it may be concluded that the value of ks3 Location 1 Location 2ShaftSpring106 satisfies a reasonable and acceptable safety limit for the deformation of the blades tip. When the stiffness was increased to more than the ks3 value as in the case of ks4, it was noted that the blade profile is far from the shroud surface and becomes closer to the base surface ks4 as seen in Fig. 5-62 (d). For the stiffness values of ks5 and ks6, the deformation of the blades was decreased to the point where that there are a few locations with a slight deformation from the base surface, as shown in Fig. 5-62 (e) and Fig. 5-62 (f), respectively. For more clarification about the stiffness variation effect, the blade profiles for all stiffness cases were combined in one graph for the most affected part in zone 2, as shown in Fig. 5-63. It may be noted from Fig. 5-63 that increasing the stiffness has a positive effect on decreasing the deformation and shifting the blade profile towards the base surface profile. It may also be noted that using a stiffness value lower than ks3 caused the blades to approach the shroud surface closely. Figure 5-62 Deformation of the blade tip surface at zone 2 for different stiffness values 107 Figure 5-63 Deformation of the blade tip surface at a specific part of zone 2 for different stiffness values 108 CHAPTER 6: CONCLUSIONS & FUTURE WORK 6.1 Conclusions 6.1.1 Conclusions on the Simulations of Rotating Stall for the Compressor With Vaned Diffuser The inception and development of rotating stall over time may be described in consecutive stages, as follows: - At a low flow rate operating condition, the flow separation occurs at the full and splitter blades’ leading edge. The random flow motion at the blades tip area causes the formation of the tip leakage flow. - The tip leakage flow causes more disturbances to the flow path inside the impeller passages and leads to an expansion of the stall area. As the stall area expands, the flow separation and tip leakage flow begin to appear at more locations, especially around the splitter blades leading edge region. The impact between the tip leakage flow and the regular inlet flow results in the formation of the tip leakage vortices. - The strength of the tip leakage vortices increases with time and these vortices move backward towards the impeller inlet. - The impeller passages become partially blocked by stall areas and the flow angle at the impeller exit deviates from the ideal value. As a result, the flow reaches the diffuser inlet with a flow angle different from the diffuser vane angle, and part of the flow fails to enter the diffuser correctly and turns back through the vaneless area; this initiates a reversed flow. It was noted that the maximum pressure fluctuations during stall occur at the blades tip area due to the effect of the tip leakage vortices; maximum pressure fluctuations also occur at the impeller exit as a result of the strong velocity variation over time at the vaneless space. The key point of the rotating stall process in centrifugal compressors is that the stall initiates at the impeller inlet due to flow separations, while stall development starts at the vaneless area due to the failure of the impeller to deliver the flow correctly to the diffuser, leading to back flow from the diffuser. 109 Stall may cause a surge as the stall area increases and the reversed flow appears at many locations at the diffuser inlet. 6.1.2 Notes on the Difference Between Compressors With Vaned and Vaneless Diffusers in the Stall Development For a compressor with a vaneless diffuser, the stall extends circumferentially at the blades tip as a thin layer that covers most of the impeller passages. This layer expands radially with time from the impeller shroud surface to the hub surface. For a compressor with a vaned diffuser, there are some impeller passages with severe stall cells and other passages that are not stalled. As time passes during stall development, stall cells start to merge and interfere with each other, forming larger stall regions. In general, at the early stages of stall, the impeller passages have a nearly identical stall effect and distribution for a compressor with a vaneless diffuser, but this distribution differs completely among passages for a compressor with a vaned diffuser, depending on which passages are affected by the back flow from the diffuser. It was found that the amplitude of pressure fluctuations at the diffuser inlet for a compressor with a vaneless diffuser is larger than that for a vaned diffuser compressor. Also, for the same time interval, the surge takes place faster in the case of a vaned diffuser compared to the case of a vaneless diffuser. 6.1.3 Comments on the Air Injection Simulations The air injection angle (Ɵ) and the injection flow rate were found to be the most powerful parameters that can determine the effectiveness of the air injection process. Compressor stability during stall can be enhanced by optimizing these parameters. It was found that the number of stall cells can be decreased by increasing the injection flow rate, except at the specific injection angles of Ɵ = 10◦ and 180◦. When the air is injected at Ɵ = 10◦, the injected air is most likely impacting with the diffuser vanes rather than being directed towards the diffuser passages. Also, the injection with an 110 angle of Ɵ = 180◦, in a direction against the flow motion, corrects the flow angle but causes flow deceleration, especially at high injection flow rates. By orienting the injected air stream closer to the impeller shroud surface by an angle of α = 5◦, the obtained results were slightly improved. The monitored pressure fluctuations at the vaneless region over time were reduced by the injection at α = 5◦, relative to the results for the angles of α = 10◦ and 15◦. This result indicates that the stall problem can be better resolved by concentrating on the regions which are very close to the shroud surface, more than other areas. The velocity distribution at the diffuser indicated that injecting air at an angle Ɵ in the range between 20◦ and 30◦ with an injection flow rate of 1.5% of the design inlet mass flow rate is the best option to correct the flow angle at the diffuser inlet and prevent the formation of reversed flow. Pressure contour results show that there are relatively high pressure areas at the impeller exit, including parts of the vaneless region, in the case of using a compressor without air injection control near the stall condition. It was found that setting the injecting angle Ɵ close to the values of 20◦ and 30◦ prevented the appearance of these high pressure regions and made the pressure distributions for most of the impeller passages identical and normal. Injection with an angle Ɵ of 30◦ provided the most stable conditions for the compressor compared to the results of other injection angles, in terms of minimizing the unsteady forces which act on the impeller blades during stall. The injection angle Ɵ of 30◦ is equivalent to nearly double of the diffuser inlet vane angle. The injection at an angle less than this value does not provide enough increase in the kinetic energy to overcome stall, while an injection at an angle larger than 30◦ value means that the injected flow starts to turn against the main flow motion, causing weakness of the main flow. For the NASA CC3 compressor, an injection at the angle Ɵ of 30◦ with an injection flow rate of 1.5% of the design inlet flow rate made the flow completely stable at an operating condition corresponding to a mass flow rate of 3.8 kg/s, and this result proves the success of the control method to shift the surge limit and increase the compressor’s stable range. 111 6.1.4 Comments on the Casing Treatment Simulations The casing grooves’ performance can be optimized by selecting the proper groove shape, size, location, and number of grooves: It was found that if the groove height is too short with respect to its width, or if the aspect ratio (height to width ratio) is less than 1, the recirculation flow generated from the groove is weak and is not sufficient to make a change in the stalled regions. If the groove cross section aspect ratio is equal to 1, the groove flow moves uniformly inside the groove and the reinjected flow velocity is enhanced compared to that for an aspect ratio of less than 1. If the groove cross section aspect ratio is larger than 1, there is no significant increase in the recirculation flow velocity compared to the case with an aspect ratio of 1, but the total area at which the flow is injected from groove increases. By locating the groove at the impeller full blades leading edge, there is a significant decrease in stall areas and blade loading. This effect diminishes gradually as the groove is located further away from this location. It was found that locating the groove at the impeller splitter blades leading edge has a relatively moderate effect on minimizing stall while locating the groove near the impeller exit has nearly no effect on minimizing stall. As the number of grooves increases, the low velocity and separation regions at the shroud surface decrease in size and move towards the impeller blades' trailing edge. On the other hand, increasing the number of grooves causes the isentropic efficiency to decrease, especially when the number of grooves exceeds 3. Groove performance at the best efficiency point is different from that during stall. At the best efficiency point, the groove flow is relatively weak compared to the main flow, and the reinjection groove flow has a nearly negligible effect on the main flow. But during stall, there is a tip leakage flow that moves against the groove flow motion, and if the tip leakage flow is stronger than the groove flow, the stall will develop till surge. If the contrary occurs and the groove flow strength is higher relatively, then the groove flow overcomes the tip leakage flow and the flow is redirected in the required direction in the impeller passages. 112 6.1.5 Conclusions on the Structural Analysis Work The amplitude of the dynamic pressure fluctuations at the blades tip increases with time as the stall develops to surge due to the interaction between the main flow and the secondary flow. The strongest fluctuations occur at the impeller tip exit due to the effect of reversed flow from the diffuser during surge. At an operating condition close to surge, the compressor is exposed to high pressure fluctuations which affect the impeller blades and cause high stresses and deformation of the blades at the tip region. The compressor model with the use of air injection provides more stability than the uncontrolled compressor model in terms of minimizing the forces and moment variation with time during stall. The most critical component which has the maximum deformation was found to be the blades tip area when the blades tip profile turns from axial to radial. This is because when the radial and axial deflection components at this area combine together, the resultant component direction is towards the shroud surface; however, for other areas such as at the impeller inlet, the blade tip profile is axial, and so the axial deformation component is not effective. It is recommended that the deflection calculations at the surge condition are taken into consideration before the final manufacturing of the compressor. It was found that the blades’ deflection during surge can cause compressor damage because the initial design for the tip clearance, designed based on steady state calculations, is smaller than the deflection during surge at some locations. By increasing the compressor shaft bearing stiffness to a specified value, it was found that the blades’ deflection at the tip area can be minimized. It would be recommended to consider this factor regarding the tip clearance value in the initial compressor design, in addition to the traditional blade deflection calculations. 113 6.2 Future Work Making more detailed optimization for the parameters of the air injection control and the casing treatment control methods. This may be infeasible because of the intensive requirement of the CPU time. According, an investigation of performing design sensitivity analysis should be considered. Regarding the air injection control, it is recommended to perform more simulations on various compressors having different geometries and then generating dimensionless relationships between the air injection parameters and the main dimensions of the compressor especially at the vanless region where the injection is recommended to be applied. Developing a design sensitivity FEM model to investigate the optimum behavior for parameter changes. Performing fully coupled fluid/structure interaction analysis and reporting the blades deflections during stall/surge. It is expected that the dynamic deformation will change the incident angle on the aerofoil and accordingly will change the dynamic behavior. Studying the compressor behavior at relatively high flow rate operating conditions close to the choke limit for a compressor enhanced by air injection control and also for a compressor modified by casing treatments. 114 REFERENCES [1] Stein, A., Niazi, S., and Sankar, L., 2000, "Computational analysis of stall and separation control in centrifugal compressors," Journal of Propulsion and Power, 16(1), pp. 65-71. [2] Greitzer, E. M., 1976, "Surge and rotating stall in axial flow compressors—Part II: experimental results and comparison with theory," Journal of Engineering for Gas Turbines and Power, 98(2), pp. 199-211. [3] Ciorciari, R., Lesser, A., Blaim, F., and Niehuis, R., 2012, "Numerical investigation of tip clearance effects in an axial transonic compressor," Journal of Thermal Science, 21(2), pp. 109-119. [4] Schleer, M., and Abhari, R. S., 2008, "Clearance effects on the evolution of the flow in the vaneless diffuser of a centrifugal compressor at part load condition," Journal of Turbomachinery, 130(3), p. 031009. [5] Tan, C., Day, I., Morris, S., and Wadia, A., 2010, "Spike-type compressor stall inception, detection, and control," Annual review of fluid mechanics, 42, pp. 275-300. [6] Haynes, J. M., Hendricks, G. J., and Epstein, A. H., 1994, "Active stabilization of rotating stall in a three-stage axial compressor," Journal of Turbomachinery, 116(2), pp. 226-239. [7] McDougall, N. M., Cumpsty, N., and Hynes, T., 1990, "Stall inception in axial compressors," Journal of Turbomachinery, 112(1), pp. 116-123. [8] Tryfonidis, M., 1994, "A study of the rotating-stall inception in high-speed compressors," Massachusetts Institute of Technology. [9] Day, I., 1993, "Stall inception in axial flow compressors," Journal of turbomachinery, 115(1), pp. 1-9. [10] Paduano, J. D., Greitzer, E., and Epstein, A., 2001, "Compression system stability and active control," Annual review of fluid mechanics, 33(1), pp. 491-517. [11] de Jager, B., "Rotating stall and surge control: A survey," Proc. Decision and Control, 1995., Proceedings of the 34th IEEE Conference on, IEEE, pp. 1857-1862. [12] Haynes, J. M., 1993, "Active control of rotating stall in a three-stage axial compressor," Massachusetts Institute of Technology. [13] Dussourd, J., Pfannebecker, G., and Singhania, S., 1977, "An experimental investigation of the control of surge in radial compressors using close coupled resistances," Journal of Fluids Engineering, 99(1), pp. 64-74. [14] Whitfield, A., Wallace, F., and Atkey, R., "The effect of variable geometry on the operating range and surge margin of a centrifugal compressor," Proc. Gas Turbine Conference and Products Show. [15] Hartmann, M. J., Benser, W. A., Hauser, C. H., and Ruggeri, R. S., 1970, "Fan and compressor technology," Paper No. NASA SP-259. [16] Osborn, W. M., Lewis Jr, G. W., and Heidelberg, L. J., 1971, "Effect of several porous casing treatments on stall limit and on overall performance of an axial flow compressor rotor." [17] Seitz, P. A., 1999, "Casing treatment for axial flow compressors," University of Cambridge. [18] Hembera, M., Kau, H.-P., and Johann, E., 2008, "Simulation of casing treatments of a transonic compressor stage," International Journal of Rotating Machinery, 2008. [19] Park, C.-Y., Choi, Y.-S., Lee, K.-Y., and Yoon, J.-Y., 2012, "Numerical study on the range enhancement of a centrifugal compressor with a ring groove system," Journal of mechanical science and technology, 26(5), pp. 1371-1378. 115 [20] Xu, W., Wang, T., and Gu, C., 2011, "Performance of a centrifugal compressor with holed casing treatment in the large flowrate condition," Science China Technological Sciences, 54(9), pp. 2483-2492. [21] Krivitzky, E., 2012, "The quest for better turbocharger compressors," Automotive Engineering International, 120(3). [22] Skoch, G. J., "Experimental investigation of centrifugal compressor stabilization techniques," Proc. ASME Turbo Expo 2003, collocated with the 2003 International Joint Power Generation Conference, American Society of Mechanical Engineers, pp. 765-776. [23] Chen, J.-P., Webster, R. S., Hathaway, M. D., Herrick, G. P., and Skoch, G. J., 2009, "High performance computing of compressor rotating stall and stall control," Integrated Computer-Aided Engineering, 16(1), pp. 75-89. [24] Hiller, S.-J., Matzgeller, R., and Horn, W., 2011, "Stability enhancement of a multistage compressor by air injection," Journal of Turbomachinery, 133(3), p. 031009. [25] Vo, H. D., Tan, C. S., and Greitzer, E. M., 2008, "Criteria for spike initiated rotating stall," Journal of turbomachinery, 130(1), p. 011023. [26] Pullan, G., Young, A., Day, I., Greitzer, E., and Spakovszky, Z., 2015, "Origins and Structure of Spike-Type Rotating Stall," Journal of Turbomachinery, 137(5), p. 051007. [27] Huang, W., Geng, S., Zhu, J., and Zhang, H., 2007, "Numerical simulation of rotating stall in a centrifugal compressor with vaned diffuser," Journal of Thermal Science, 16(2), pp. 115-120. [28] Geng, S., Lin, F., Chen, J., and Nie, C., 2011, "Evolution of unsteady flow near rotor tip during stall inception," Journal of Thermal Science, 20(4), pp. 294-303. [29] Yamada, K., Kikuta, H., Iwakiri, K.-i., Furukawa, M., and Gunjishima, S., 2013, "An explanation for flow features of spike-type stall inception in an axial compressor rotor," Journal of turbomachinery, 135(2), p. 021023. [30] Tomita, I., Ibaraki, S., Furukawa, M., and Yamada, K., 2013, "The Effect of Tip Leakage Vortex for Operating Range Enhancement of Centrifugal Compressor," Journal of Turbomachinery, 135(5), p. 051020. [31] Ramakrishna, P., and Govardhan, M., 2009, "Stall characteristics and tip clearance effects in forward swept axial compressor rotors," Journal of Thermal Science, 18(1), pp. 40-47. [32] Vagani, M., Engeda, A., and Cave, M., 2013, "Prediction of impeller rotating stall onset using numerical simulations of a centrifugal compressor. Part 1: Detection of rotating stall using fixed-flow transient simulations," Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, p. 0957650912474386. [33] Bailey, E. E., 1972, "Effect of grooved casing treatment on the flow range capability of a single-stage axial-flow compressor." [34] Huang, X., Chen, H., and Fu, S., "CFD investigation on the circumferential grooves casing treatment of transonic compressor," Proc. XIX International Symposium on Air Breathing Engines (ISABE 2009), Montreal, Canada, September, pp. 7-11. [35] Voges, M., Schnell, R., Willert, C., Mönig, R., Müller, M., and Zscherp, C., 2011, "Investigation of blade tip interaction with casing treatment in a transonic compressor—Part I: particle image velocimetry," Journal of Turbomachinery, 133(1), p. 011007. [36] Kurokawa, J., Saha, S. L., Matsui, J., and Kitahora, T., 2000, "Passive control of rotating stall in a parallel-wall vaneless diffuser by radial grooves," Journal of fluids engineering, 122(1), pp. 90-96. 116 [37] Kim, J.-H., Choi, K.-J., and Kim, K.-Y., 2013, "Aerodynamic analysis and optimization of a transonic axial compressor with casing grooves to improve operating stability," Aerospace Science and Technology, 29(1), pp. 81-91. [38] Xi, G., Ma, Y., Wu, G., Zhang, K., Xiao, W., and Mou, Z., 2013, "Exploration of a new-type grooved casing treatment configuration for a high-speed small-size centrifugal compressor," Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, p. 0957650913477094. [39] Houghton, T., and Day, I., 2012, "Stability enhancement by casing grooves: The importance of stall inception mechanism and solidity," Journal of Turbomachinery, 134(2), p. 021003. [40] Gelmedov, F. S., Lokshtanov, E. A., Olstain, L. E.-M., and Sidorkin, M. A., 1998, "Anti-stall tip treatment means," Google Patents. [41] Prince, D., Wisler, D., and Hilvers, D., "A study of casing treatment stall margin improvement phenomena," Proc. ASME 1975 International Gas Turbine Conference and Products Show, American Society of Mechanical Engineers, pp. V01AT01A059-V001AT001A059. [42] Lu, X., Chu, W., Zhu, J., and Zhang, Y., 2009, "Numerical investigations of the coupled flow through a subsonic compressor rotor and axial skewed slot," Journal of Turbomachinery, 131(1), p. 011001. [43] Xu, W., Wang, T., Gu, C., and Ding, L., 2012, "Numerical Investigation of a Centrifugal Compressor With Holed Casing Treatment," Journal of engineering for gas turbines and power, 134(4), p. 044502. [44] Ding, L., Wang, T., Yang, B., Xu, W., and Gu, C., 2013, "Experimental investigation of the casing treatment effects on steady and transient characteristics in an industrial centrifugal compressor," Experimental Thermal and Fluid Science, 45, pp. 136-145. [45] Yang, M., Martinez-Botas, R., Zhang, Y., Zheng, X., Bamba, T., Tamaki, H., and Li, Z., "Investigation of self-recycling-casing-treatment (srct) influence on stability of high pressure ratio centrifugal compressor with a volute," Proc. ASME 2011 Turbo Expo: Turbine Technical Conference and Exposition, American Society of Mechanical Engineers, pp. 1857-1867. [46] Skoch, G. J., "Experimental investigation of diffuser hub injection to improve centrifugal compressor stability," Proc. ASME Turbo Expo 2004: Power for Land, Sea, and Air, American Society of Mechanical Engineers, pp. 793-804. [47] Lu, X., Chu, W., Zhu, J., and Tong, Z., 2006, "Numerical and Experimental Investigations of Steady Micro-Tip Injection on a Subsonic Axial-Flow Compressor Rotor," International Journal of Rotating Machinery, 2006. [48] Khaleghi, H., Boroomand, M., Teixeira, J., and Tousi, A., 2008, "A numerical study of the effects of injection velocity on stability improvement in high-speed compressors," Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 222(2), pp. 189-198. [49] Sadeghi, M., and Liu, F., "Coupled fluid-structure simulation for turbomachinery blade rows," Proc. 43rd AIAA Aerospace Sciences Meeting and Exhibit, pp. 2005-0018. [50] Clark, W. S., and Hall, K. C., "A time-linearized Navier-Stokes analysis of stall flutter," Proc. ASME 1999 International Gas Turbine and Aeroengine Congress and Exhibition, American Society of Mechanical Engineers, pp. V004T003A040-V004T003A040. [51] Carstens, V., and Belz, J., 2001, "Numerical investigation of nonlinear fluid-structure interaction in vibrating compressor blades," Journal of turbomachinery, 123(2), pp. 402-408. 117 [52] Johnston, D. A., Cross, C. J., and Wolff, J. M., 2005, "An architecture for fluid/structure interaction analysis of turbomachinery blading," AIAA paper no, 4013, pp. 10-13. [53] Gnesin, V., Rzadkowski, R., and Kolodyazhnaya, L., 2000, "A Coupled Fluid-Structure Analysis for 3D Flutter in Turbomachines," ASME J. [54] Srivastava, R., Bakhle, M., Keith Jr, T., and Stefko, G., 2001, "Aeroelastic stability computations for turbomachinery." [55] Im, H.-S., and Zha, G.-C., "Simulation of non-synchronous blade vibration of an axial compressor using a fully coupled fluid/structure interaction," Proc. ASME Turbo Expo 2012: Turbine Technical Conference and Exposition, American Society of Mechanical Engineers, pp. 1395-1407. [56] Brandsen, J. D., 2013, "Prediction of axial compressor blade vibration by modelling fluid-structure interaction," Stellenbosch: Stellenbosch University. [57] Sadeghi, M., and Liu, F., 2005, "Computation of cascade flutter by uncoupled and coupled methods," International Journal of Computational Fluid Dynamics, 19(8), pp. 559-569. [58] Pei, J., Yuan, S., and Yuan, J., 2013, "Fluid-structure coupling effects on periodically transient flow of a single-blade sewage centrifugal pump," Journal of Mechanical Science and Technology, 27(7), pp. 2015-2023. [59] FLUENT, A., 2011, "13.0. ANSYS," Inc., Canonsburg, PA, USA. [60] Fluent, A., 2011, "Ansys Fluent Theory Guide," ANSYS Inc., USA. [61] Davis, P., Rinehimer, A., and Uddin, M., "A comparison of RANS-based turbulence modeling for flow over a wall-mounted square cylinder," Proc. 20th Annual Conference of the CFD Society of Canada. [62] Wilcox, D. C., 1998, Turbulence modeling for CFD, DCW industries La Canada, CA. [63] Hinze, J., 1975, "TurbulenceMcGraw-Hill," New York, 218. [64] Pope, S. B., 2000, Turbulent flows, Cambridge university press. [65] Launder, B. E., and Spalding, D., 1974, "The numerical computation of turbulent flows," Computer methods in applied mechanics and engineering, 3(2), pp. 269-289. [66] Jayatilleke, C. L. V., 1966, "The influence of Prandtl number and surface roughness on the resistance of the laminar sub-layer to momentum and heat transfer," University of London. [67] Li, R., Tang, T., and Zhang, P., 2002, "A moving mesh finite element algorithm for singular problems in two and three space dimensions," Journal of Computational Physics, 177(2), pp. 365-393. [68] Murayama, M., Nakahashi, K., and Matsushima, K., 2002, "Unstructured dynamic mesh for large movement and deformation," AIAA paper, 122, p. 2002. [69] Skoch, G. J., Prahst, P., Wernet, M., Wood, J., and Strazisar, A., "Laser anemometer measurements of the flow field in a 4: 1 pressure ratio centrifugal impeller," Proc. ASME 1997 International Gas Turbine and Aeroengine Congress and Exhibition, American Society of Mechanical Engineers, pp. V001T003A049-V001T003A049. [70] Spakovsky, Z., "Backward traveling rotating stall waves in centrifugal compressors," Proc. ASME Turbo Expo 2002: Power for Land, Sea, and Air, American Society of Mechanical Engineers, pp. 529-543. [71] Tarr, D. L., 2008, "Scaling of impeller response to impeller-diffuser interactions in centrifugal compressors," Massachusetts Institute of Technology. [72] Wernet, M. P., Bright, M. M., and Skoch, G. J., 2001, "An investigation of surge in a high-speed centrifugal compressor using digital PIV," Journal of turbomachinery, 123(2), pp. 418-428. 118 [73] McKain, T. F., 1997, Coordinates for a high performance 4: 1 pressure ratio centrifugal compressor, National Aeronautics and Space Administration. [74] Skoch, G. J., 2000, "Centrifugal Compressor Flow Range Extension Using Diffuser Flow Control," NASA Glenn Research Center, Cleveland, Ohio, Dec, 5, p. 17. [75] ANSYS, R., 2011, "14.0, Help System,“Coupled Field Analysis Guide”, ANSYS," Inc. [76] Richardson, L. F., 1911, "The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam," Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, pp. 307-357. [77] Larosiliere, L., Skoch, G. J., and Prahst, P., 1997, Aerodynamic synthesis of a centrifugal impeller using CFD and measurements, National Aeronautics and Space Administration. [78] Nie, C., Tong, Z., Geng, S., Zhu, J., and Huang, W., 2007, "Experimental investigations of micro air injection to control rotating stall," Journal of Thermal Science, 16(1), pp. 1-6. [79] Spakovszky, Z., and Roduner, C., 2009, "Spike and modal stall inception in an advanced turbocharger centrifugal compressor," Journal of Turbomachinery, 131(3), p. 031012. [80] While, M., 1979, "Rolling element bearing vibration transfer characteristics: effect of stiffness," Journal of applied mechanics, 46(3), pp. 677-684. [81] Guo, Y., and Parker, R. G., 2012, "Stiffness matrix calculation of rolling element bearings using a finite element/contact mechanics model," Mechanism and Machine Theory, 51, pp. 32-45. 119 APPENDIX A.1 The Conservation Equation of Transport and its Discretization The discretization can be described by considering the conservation equation of transport for the scalar variable ∅ for an arbitrary control volume ܸ as shown in equation A-1 which is applied to each control volume. The discretization of equation A-1 on a given cell results in equation A-2. රρϕv.ሬሬԦ ݀ܣሬሬሬሬሬԦ ൌ රΓம∅. ݀ܣሬሬሬሬሬԦ න ܵ∅ܸ݀ (A-1) Where ߩ is the fluid density, vሬԦ is the velocity vector (vሬԦ ൌ ݑଓ̂ ݒଔ̂ ݓ ݇ ), ܣԦis the surface area vector,Γ∅ is the diffusion coefficient for ∅ and ܵ∅ is the source of ∅ per unit volume. Also, ∅ൌ ߲∅߲ݔ ଓ̂ ߲∅߲ݕ ଔ̂ ߲∅߲ݖ ݇. ߩݒሬሬሬሬԦ∅ ൌ Γமሺ∅ሻ. ܣሬሬሬሬԦ ܵ∅ܸ (A-2) Where ݊ is the number of faces enclosing the cell, ∅ is the value of ∅ convected through the face f, ߩݒሬሬሬሬԦ∅ is the mass flux through the face, ܣሬሬሬሬԦ is the area of face f for (|ܣ| ൌ หܣ௫ଓ̂ ܣ௬ଔ̂ ܣ௭ ݇ห), ሺ∅ሻ is the magnitude of ∅ normal to face f, and ܸ is the cell volume. A.2 Boussinesq Hypothesis Equation A-3 shows the concept of the Boussinesq hypothesis. െߩݑప́ ݑఫ́തതതതത ൌ ߤ௧ ቆ߲ݑ߲ݔ ߲ݑ߲ݔቇ െ23 ሺߩ݇ ߤ௧߲ݑ߲ݔሻߜ (A-3) Where ߤ௧ is the turbulent viscosity which is computed by solving the transport equations for the turbulence kinetic energy (݇) and the turbulence dissipation rate (ߝ). 120 A.3 Realizable െ ࢿ Model Equations Equation A-4 indicates the method of calculation of the turbulent viscosity ߤ௧ for the realizable ݇ െ ߝ turbulence model. ߤ௧ ൌ ߩܥఓ ݇ଶߝ (A-4) Where ܥఓis calculated as shown in equation 3-5. ܥఓ ൌ 14.04 ఌ ܤ (A-5) Where ܤ is a parameter which depends on many factors related to velocity gradients and shear stresses as described in equation A-6. ܤ ൌ ሾට6൫ ܵ ܵ Ωపఫതതതത൯ሿሾcosሺ13 cosିଵሺ√6 ܵ ܵܵܵ̅ଷ ሻሻሿ (A-6) Where ܵ is the shear stress and ܵ̅ refers to the average shear stress while Ωഥ represents the mean rate of rotation tensor. Equations A-7 and A-8 indicate how ܵ̅ and ܵ are calculated respectively. ܵ̅ ൌ ට ܵ ܵ (A-7) ܵ ൌ 12 ቆ߲ݑ߲ݔ ߲ݑ߲ݔቇ (A-8) A.4 Definitions of the Near-Wall Treatment Parameters The dimensionless wall distance ݕା is defined as shown in equation A-9. ݕା ൌ ߩݕ ఛܷμ (A-9) Where ݕ is the distance from point P to the wall, ߤ is the dynamic viscosity of the fluid, ߩ is the fluid density, while the definition of ఛܷ (the frictional velocity) is given in equation A-10. ఛܷ ൌ ඨ߬௪ߩ (A-10) Where ߬௪ is the wall shear stress and it is computed as shown in equation A-11. 121 ߬௪ ൌ μ߲ܷ߲ݕ (A-11) Where ܷ is the mean velocity of the fluid at point P Also, the ܷା is the dimensionless velocity and its definition is as presented in equation A-12. ܷା ൌ ܷఛܷ (A-12) The logarithmic law of the wall is applied at the inner layer region and it is described in equation A-13. ܷା ൌ 1ߢ ln ݕା Cା (A-13) Where Cା is a constant depends on the wall surface roughness and ߢ is the von Kármán constant (= 0.42). The logarithmic law of the wall can be used for ݕା 11 When ݕା ൏ 11, then ܷା is nearly equal to ݕା. A.5 Equations for the Using of Wall Functions for Momentum The dimensionless parameters ܷ∗ and ݕ∗ are described in equations A-14 & A-15 respectively. ܷ∗ ൌ ܷܥఓଵ ସ⁄ ݇ଵ ଶ⁄߬௪ ߩ⁄ (A-14) ݕ∗ ൌ ߩܥఓଵ ସ⁄ ݇ଵ ଶ⁄ ݕߤ (A-15) Where ܷ is the mean velocity of the fluid at point P, ݇ is the turbulence kinetic energy at point P, ݕ is the distance from point P to the wall, and ߤ is the dynamic viscosity of the fluid. The law of the wall for the mean velocity can be described as shown in equation A-16. (A-16) Where ߢ is the von Kármán constant (= 0.42), ܧ is an empirical constant (= 9.793). 122 When ݕ∗< 11.225 at the cells adjacent to wall, the laminar relation is used as shown in equation A-17. (A-17) A.6 Equations for the Using of Wall Functions for Energy The law of the wall for temperature is described in equation A-18. ܶ∗ ൌ ሺ ௪ܶ െ ܶሻߩܿܥఓଵ ସ⁄ ݇ଵ ଶ⁄ݍሶ ܶ∗ ൌۖە۔ۖۓ Pr ݕ∗ 12ߩ Prܥఓଵ ସ⁄ ݇ଵ ଶ⁄ ܷଶݍሶ ሺݕ∗ ൏ ݕ∗் ሻܲݎ௧ 1ߢ lnሺܧݕ∗ሻ ܲ൨ 12ߩPrܥఓଵ ସ⁄ ݇ଵ ଶ⁄ݍሶ ሼܲݎ௧ܷଶ ሺPr െ ܲݎ௧ሻ ܷଶሽ ሺݕ∗ ݕ∗்ሻ (A-18) Where P is calculated based on equation A-19 (A-19) Where ܿ represents the specific heat of the fluid, ݍሶ is the wall heat flux, ܶ describes the temperature at the cell adjacent to wall, ௪ܶ is the wall temperature, ܲݎ is the molecular Prandtl number, ܲݎ௧ is the turbulent Prandtl number and ܷ is the mean velocity magnitude when ݕ∗ ൌݕ∗் . The thermal sub-layer thickness ݕ∗் presented in equation A-18 is calculated like ݕ∗ value for which the logarithmic and linear laws intersect resulting in the definition of the fluid molecular Prandtl number. The application of the law of the wall for temperature has a specific procedure. Based on the values of the fluid physical properties, the molecular Prandtl number is calculated. After that, the thermal sub layer thickness ݕ∗் is computed by intersecting the linear and logarithmic profiles. Based on the ݕ∗ value near wall, the linear or the logarithmic profile is applied in order to compute the wall temperature ௪ܶ or the heat flux ݍሶ which depends on the thermal boundary condition type.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical investigation of rotating stall in centrifugal...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical investigation of rotating stall in centrifugal compressors Halawa, Taher Mohammed Abdelhakim Fattouh 2016
pdf
Page Metadata
Item Metadata
Title | Numerical investigation of rotating stall in centrifugal compressors |
Creator |
Halawa, Taher Mohammed Abdelhakim Fattouh |
Publisher | University of British Columbia |
Date Issued | 2016 |
Description | Rotating stall is a complicated phenomenon which may occur when a centrifugal compressor works at low flow rate condition close to the compressor surge limit. Rotating stall and surge may cause a large amount of damage to the compressor components due to instability and high pressure fluctuations which increase vibration of the compressor blades. This research introduces numerical simulations of the rotating stall phenomenon in a high speed centrifugal compressor with vaned and vaneless diffusers. The present study discusses the improvement of two important methods to control stall: the air injection method and the casing grooves treatment method. The method of air injection was optimized by monitoring the effects resulting from the variation of the angle of injection and the injection flow rate. The efficiency of the casing grooves method was related to the grooves' locations, cross section dimensions, and the number of grooves. The current research also presents estimations of the deformation of the compressor blades during the transition stage between stall and surge in order to determine whether the blades will deform safely within the tip clearance gap without touching the casing. Rotating stall simulations showed that the initiation of stall is due to the flow separation and the formation of the tip leakage flow. Simulations showed also that stall development appears at the vaneless area between the impeller and diffuser. Air injection simulations showed that stability is optimized when the air injection angle is approximately double the value of the diffuser inlet vane angle. Simulations showed that in order to achieve effective flow circulation, the aspect ratio of the casing groove cross section should be at least equal to one. Structural analysis showed that the most dangerous region at the deformed blade’s tip surface is located at the transition area between the axial and radial blade profile. Results also showed that the deformation of the blades at the tip region can be limited by changing the shaft bearing stiffness; this is useful to prevent damage to the compressor blades. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2016-09-22 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0314571 |
URI | http://hdl.handle.net/2429/59261 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2016-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2016_november_halawa_taher.pdf [ 16.99MB ]
- Metadata
- JSON: 24-1.0314571.json
- JSON-LD: 24-1.0314571-ld.json
- RDF/XML (Pretty): 24-1.0314571-rdf.xml
- RDF/JSON: 24-1.0314571-rdf.json
- Turtle: 24-1.0314571-turtle.txt
- N-Triples: 24-1.0314571-rdf-ntriples.txt
- Original Record: 24-1.0314571-source.json
- Full Text
- 24-1.0314571-fulltext.txt
- Citation
- 24-1.0314571.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:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0314571/manifest