INVESTIGATION OF THE PERFORMANCE OF LABYRINTH SEALS FOR CENTRIFUGAL COMPRESSOR APPLICATIONS By Casey Palanca A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering—Doctor of Philosophy 2022 ABSTRACT INVESTIGATION OF THE PERFORMANCE OF LABYRINTH SEALS FOR CENTRIFUGAL COMPRESSOR APPLICATIONS By Casey Palanca Labyrinth seals were one of the first seal configurations used in modern turbomachinery and, due to their robust qualities and relatively low-cost productions, remain one of the most widely used seal configurations today. Their primary purpose is to control internal leakage between the rotating and stationary components of rotating machinery, including the centrifugal compressor. A reduction in secondary leakage flow will always be accompanied by an increase in efficiency. However, while fulfilling the objective of restricting secondary fluid flow, labyrinth seals have been known to cause adverse stability effects on the rotor. Driving forces inside the cavities from the circumferential flow path have been known to be a potential source of destabilizing vibrations. Therefore, accurately predicting these forces is a primary interest in compressor design. These forces are characterized by stiffness and damping coefficients. The present study utilizes the growing advances in CFD to understand, model, and predict the aerodynamic and rotordynamic performance of labyrinth seals. The scope of this work progresses from a well-established steady state CFD method to a more novel transient CFD approach. The benefits and disadvantages of each method are evaluated and discussed by comparing accuracy, reliability, and computational efficiency. Each method is validated with experimental data. Additionally, the proposed transient CFD method can be used to perform a reasonably accurate prediction of the frequency-dependent rotordynamic coefficients by using a Fast Fourier Transform analysis on the monitored force response and displacement data. Lastly, the transient CFD approach is expanded upon by investigating the flow characteristics of long 18 tooth on rotor balance piston labyrinth seal modeled with abradable grooves on the stator. It was discovered that the creation of abradable grooves on the stator can cause the vortex between the labyrinth teeth to change directions (clockwise to counterclockwise). This observation is used to determine the relationship between the flow pattern and rotordynamic performance. A parametric study shows the effect of abradable groove geometries and operating flow conditions on the labyrinth seal rotordynamic coefficients. ACKNOWLEDGMENTS There are several people that I need to acknowledge without whom this work would not have been possible. First, my Lord and Savior Jesus Christ. “For from him and through him and for him are all things. To him be the glory forever! Amen”. Second, I would like to express tremendous gratitude to my advisor Dr. Abraham Engeda for choosing me to work on this project and supporting me along the way. I couldn’t have asked for a better advisor. I would also like to thank Dr. Norbert Mueller, Dr. Andre Benard, and Dr. Wei Liao for agreeing to participate in my guidance committee. This work would not be possible without the support provided by Solar Turbines Incorporated. I would like to thank Mike Cave and the aerodynamic and rotordynamic group for their support and guidance during my internships and research work. A special thanks to Marco Vagani and Chris Clarke, who helped steer me in the right direction early on in my graduate career. On a personal note, I would like to thank my family who has always been a tremendous beacon of support, and my church family and friends who have helped me preserve as I navigated each milestone iv TABLE OF CONTENTS LIST OF TABLES ..................................................................................................................... vii LIST OF FIGURES .................................................................................................................. viii KEY TO SYMBOLS ................................................................................................................ xiv Chapter 1: Background and Motivation .............................................................................. 1 1.1. Centrifugal Compressors .................................................................................................. 1 1.2. Industrial Compressors ..................................................................................................... 2 1.3. Introduction to Compressor Seals .................................................................................... 2 1.3.1. Shaft-End Seals ......................................................................................................... 3 1.3.2. Internal Labyrinth Seals ............................................................................................ 4 1.4. Gas Labyrinth Seal Applications...................................................................................... 5 1.5. Objectives ......................................................................................................................... 6 Chapter 2: Evolution of Theory and Flow Through Labyrinth Seals .................................. 9 2.1. Flow Structure and Features in the Labyrinth Seals ........................................................ 9 2.2. Rotordynamics ................................................................................................................... 12 2.2.1. Rotordynamic Coefficients ..................................................................................... 15 2.3. Literature Review ........................................................................................................... 19 Chapter 3: Methodology ......................................................................................................... 27 3.1. Experimental Benchmark ............................................................................................... 27 3.2. Computational Fluid Dynamics ..................................................................................... 29 3.3. Computational Models ................................................................................................... 30 3.3.1. Tooth on Rotor Labyrinth Seal ............................................................................... 30 3.3.2. Tooth on Stator Labyrinth Seal ............................................................................... 33 3.4. Turbulence Models ......................................................................................................... 34 3.5. Grid Generation .............................................................................................................. 38 3.6. Steady State CFD ........................................................................................................... 42 3.7. Transient CFD ................................................................................................................ 44 3.7.1. Rotor Whirling Models ........................................................................................... 46 3.7.2. Transient Rotordynamic Coefficient Identification ................................................ 50 Chapter 4: Steady State CFD Results .................................................................................... 54 4.1. 8-Tooth on Rotor Labyrinth Seal ................................................................................... 54 4.2. 5-Stepped Labyrinth Seal ............................................................................................... 58 4.3 5-Tooth on Stator Labyrinth Seal ........................................................................................ 59 4.4. Conclusions .................................................................................................................... 62 Chapter 5: Steady State vs Transient CFD Comparison ..................................................... 63 5.1. CFD Model and Set up ................................................................................................... 63 5.2. Post Processing............................................................................................................... 65 5.2.1 Computational time................................................................................................. 87 5.3. Comparison to Experiment............................................................................................. 88 5.3.1 Stability Characteristics ........................................................................................ 100 v 5.4. Conclusions .................................................................................................................. 102 Chapter 6: Transient CFD: Further Investigation ............................................................. 104 6.2. Flow Characteristics ..................................................................................................... 107 6.4. Rotordynamic Coefficients .......................................................................................... 119 6.4.1 Abradable grooves ................................................................................................ 128 6.5. Conclusions .................................................................................................................. 141 Chapter 7: Conclusions ......................................................................................................... 143 7.1. Summary and Recommendations ................................................................................. 143 REFERENCES ........................................................................................................................ 143 vi LIST OF TABLES Table 1: Boundary conditions for TOR labyrinth seals ................................................................ 31 Table 2: Boundary conditions, balance piston labyrinth seal ....................................................... 33 Table 3: Boundary conditions for TOS labyrinth seals ................................................................ 34 Table 4: Number of nodes in Figure 15 for each seal ................................................................... 40 Table 5: Rotordynamic coefficients For TOR labyrinth seal at 0 inlet swirl ratio ....................... 57 Table 6: Rotordynamic coefficients For TOR labyrinth seal at 0.165 inlet swirl ratio ................ 57 Table 7: Comparison of predicted seal coefficients with the results reported by [37, 41, 42] ..... 61 Table 8: Transient KCM CFD model computational time ........................................................... 88 Table 9: Balance Piston Dimensions .......................................................................................... 105 Table 10: Abradable groove parametric study ............................................................................ 106 Table 11: Leakage comparison, groove, with and without abradable grooves ........................... 114 Table 12: Rotor whirling frequency ............................................................................................ 119 vii LIST OF FIGURES Figure 1: Labyrinth seal configurations a) TOS b) TOR c) Interlocking ....................................... 4 Figure 2: Typical labyrinth seal leakage for a centrifugal compressor [9] ..................................... 5 Figure 3: Example of the flow path through labyrinth seal ............................................................ 9 Figure 4: shroud face configurations of impeller eye seal a) open face, b) straight face c) curved face. [3].......................................................................................................................... 12 Figure 5: The Jeffcott rotor [18] ................................................................................................... 14 Figure 6: Whirling Rotor Model [16] ........................................................................................... 14 Figure 7: Rotordynamic model of gas labyrinth seal with a centered rotor.................................. 15 Figure 8: Direct and cross-coupled seal reaction forces acting on a whirling rotor ..................... 16 Figure 9: Conventional balance piston labyrinth seal flow [21]. .................................................. 18 Figure 10: Balance piston labyrinth seal with shunt injection [21]. ............................................. 18 Figure 11: Comparison of predicted pressure distribution with experimental measurements ...... 31 Figure 12: Operating flow points along the compressor map ....................................................... 32 Figure 13: Boundary Conditions Used for the Balance Piston Simulations ................................. 33 Figure 14: Law of the Wall ........................................................................................................... 36 Figure 15: Regions of node refinement (a) TOS and TOR (b) Stepped labyrinth seal ................ 40 Figure 16: Circumferential node study for 8 TOR and 5-Stepped labyrinth seal ......................... 41 Figure 17: Circumferential node study for 5 TOS labyrinth seal ................................................. 41 Figure 18: Cavity mesh a) 8-TOR b) 5-Stepped c) 5-TOS ........................................................... 42 Figure 19: Stationary to rotating frame transfer ........................................................................... 44 Figure 20: Derivation of the Rotor Surface Velocity.................................................................... 46 Figure 21: Circular whirl orbit a) forward whirl, b) backward whirl ........................................... 48 Figure 22: Elliptical whirl orbit a) x-direction excitation, b) y-direction excitation .................... 49 Figure 23: Radial impedance for the TOR labyrinth at 0 inlet swirl ............................................ 55 Figure 24: Tangential impedance for the TOR labyrinth at 0 inlet swirl...................................... 56 Figure 25: Radial impedance for the TOR labyrinth at 0.165 inlet swirl ..................................... 56 Figure 26: Tangential impedance for the TOR labyrinth at 0.165 inlet swirl............................... 57 viii Figure 27: Comparison of rotordynamic coefficients for 5-tooth stepped seal (a) direct stiffness (b) cross-coupled stiffness (c) direct damping (d) cross-coupled damping .................. 59 Figure 28: Radial impedances for 5-TOS labyrinth at 0 inlet swirl .............................................. 61 Figure 29:Tangential impedances for 5-TOS labyrinth at 0 inlet swirl ........................................ 61 Figure 30: Comparison of Experimental Leakage measurements with steady and transient CFD ....................................................................................................................................... 64 Figure 31: Swirl ratio through each seal cavity ............................................................................ 65 Figure 32: Monitored whirling motion for circular orbit, forward excitation in the time domain 68 Figure 33: Monitored whirling motion for circular orbit, backward excitation in the time domain ....................................................................................................................................... 69 Figure 34:Monitored whirling motion for elliptical orbit, x excitation in the time domain ......... 69 Figure 35: Monitored whirling motion for elliptical orbit, y excitation in the time domain ........ 70 Figure 36: FFT of whirling motion for circular orbit, forward and backward excitation ............. 70 Figure 37: FFT of whirling motion, elliptical orbit, x excitation................................................. 71 Figure 38: FFT of whirling motion, elliptical orbit, y excitation.................................................. 71 Figure 39: Monitored force response for single-frequency circular orbit, backward excitation. 0.08 preswirl. ................................................................................................................. 72 Figure 40: Monitored force response for single-frequency circular orbit, forward excitation 0.08 preswirl. ......................................................................................................................... 72 Figure 41: Monitored force response for multi-frequency circular orbit, backward excitation. 0.08 preswirl. ................................................................................................................. 73 Figure 42: FFT of fluid force response for circular orbit, backward excitation. 0.08 preswirl ... 73 Figure 43: Monitored force response for multi-frequency circular orbit, forward excitation. 0.08 preswirl. ......................................................................................................................... 74 Figure 44: FFT of fluid force response for circular orbit, forward excitation. 0.08 preswirl ....... 74 Figure 45: Monitored force response for multi-frequency elliptical orbit, x excitation. .............. 75 Figure 46: FFT of fluid force response for elliptical, x excitation. 0.08 preswirl ........................ 75 Figure 47: Monitored force response for multi-frequency elliptical orbit, y excitation ............... 76 Figure 48: FFT of fluid force response for elliptical, y excitation. 0.08 preswirl ........................ 76 Figure 49: Monitored force response for single-frequency circular orbit, backward excitation. 0.25 preswirl. ................................................................................................................. 77 Figure 50: Monitored force response for single-frequency circular orbit, forward excitation. 0.25 preswirl. ......................................................................................................................... 77 ix Figure 51: Monitored force response for multi-frequency circular orbit, backward excitation. 0.25 preswirl. ................................................................................................................. 78 Figure 52: FFT of fluid force response for circular orbit, backward excitation. 0.25 preswirl .... 78 Figure 53: Monitored force response for multi-frequency circular orbit, forward excitation. 0.25 preswirl. ......................................................................................................................... 79 Figure 54: FFT of fluid force response for circular orbit, forward excitation. 0.25 preswirl ....... 79 Figure 55: Monitored force response for multi-frequency elliptical orbit, x excitation. 0.25 preswirl. ......................................................................................................................... 80 Figure 56: FFT of fluid force response for elliptical orbit, x excitation. 0.25 preswirl ................ 80 Figure 57: Monitored force response for multi-frequency elliptical orbit, y excitation. 0.25 preswirl. ......................................................................................................................... 81 Figure 58: FFT of fluid force response for elliptical orbit, y excitation. 0.25 preswirl ................ 81 Figure 59: Monitored force response for single-frequency circular orbit, backward excitation. 0.48 preswirl. ................................................................................................................. 82 Figure 60: Monitored force response for single-frequency circular orbit, forward excitation. 0.48 preswirl .......................................................................................................................... 82 Figure 61: Monitored force response for multi-frequency circular orbit, backward excitation. 0.48 preswirl. ................................................................................................................. 83 Figure 62: FFT of fluid force response for circular orbit, backward excitation. 0.48 preswirl .... 83 Figure 63: Monitored force response for multi-frequency circular orbit, forward excitation. 0.48 preswirl. ......................................................................................................................... 84 Figure 64: FFT of fluid force response for circular orbit, forward excitation. 0.48 preswirl ....... 84 Figure 65: Monitored force response for multi-frequency elliptical orbit, x excitation. 0.48 preswirl. ......................................................................................................................... 85 Figure 66: FFT of fluid force response for elliptical orbit, x excitation. 0.48 preswirl ................ 85 Figure 67: Monitored force response for multi-frequency elliptical orbit, y excitation. 0.48 preswirl. ......................................................................................................................... 86 Figure 68: FFT of fluid force response for elliptical orbit, y excitation. 0.48 preswirl ................ 86 Figure 69: Comparison of direct stiffness for steady state and transient KCM with experiment. 92 Figure 70: Comparison of predicted direct stiffness using two multifrequency whirling models with experiment. 0.08 preswirl. ..................................................................................... 92 Figure 71: Comparison of predicted direct stiffness using two multifrequency whirling models with experiment. 0.25 preswirl. ..................................................................................... 93 x Figure 72: Comparison of predicted direct stiffness using two multifrequency whirling models with experiment. 0.48 preswirl. ..................................................................................... 93 Figure 73: Comparison of cross coupled stiffness for steady state and transient KCM model with experiment ..................................................................................................................... 94 Figure 74: Comparison of predicted cross coupled stiffness using two multifrequency whirling models with experiment. 0.08 preswirl. ........................................................................ 94 Figure 75: Comparison of predicted cross coupled stiffness using two multifrequency whirling models with experiment. 0.25 preswirl. ........................................................................ 95 Figure 76: Comparison of predicted cross coupled stiffness using two multifrequency whirling models with experiment. 0.48 preswirl. ........................................................................ 95 Figure 77: Comparison of direct damping for steady state and transient KCM model with experiment ..................................................................................................................... 96 Figure 78: Comparison of predicted direct damping using two multifrequency whirling models with experiment. 0.08 preswirl. ..................................................................................... 96 Figure 79: Comparison of predicted direct damping using two multifrequency whirling models with experiment. 0.25 preswirl. ..................................................................................... 97 Figure 80: Comparison of predicted direct damping using two multifrequency whirling models with experiment. 0.48 preswirl. ..................................................................................... 97 Figure 81: Comparison of cross coupled damping for steady state and transient KCM model with experiment ..................................................................................................................... 98 Figure 82: Comparison of predicted cross coupled damping using two multifrequency whirling models with experiment. 0.08 preswirl. ........................................................................ 98 Figure 83: Comparison of predicted cross coupled damping using two multifrequency whirling models with experiment. 0.25 preswirl. ........................................................................ 99 Figure 84: Comparison of predicted cross coupled damping using two multifrequency whirling models with experiment. 0.48 preswirl. ........................................................................ 99 Figure 85: Comparison of predicted effective damping using two multifrequency whirling models with experiment. 0.08 preswirl. ...................................................................... 101 Figure 86: Comparison of predicted effective damping using two multifrequency whirling models with experiment. 0.25 preswirl. ...................................................................... 101 Figure 87: Comparison of predicted effective damping using two multifrequency whirling models with experiment. 0.48 preswirl. ...................................................................... 102 Figure 88: Balance Piston Labyrinth Seal Dimensions ............................................................. 105 Figure 89: Detailed Labyrinth Seal Tooth Dimensions .............................................................. 105 Figure 90: Balance piston seal mesh, section cut view ............................................................... 107 xi Figure 91: Circumferential node study for 18 TOR balance piston seal .................................... 107 Figure 92: Pressure distribution through the labyrinth seal. a) choke b) midpoint c) surge ....... 108 Figure 93: Detailed example of traditional flow pattern through labyrinth seal......................... 109 Figure 94: Velocity vectors in stationary frame ......................................................................... 110 Figure 95: Recirculation zones for TOR seal with abradable grooves ....................................... 112 Figure 96: Flow pattern for TOR seal with smooth stator .......................................................... 112 Figure 97: Circumferential velocity in the stationary frame, no abradable grooves .................. 112 Figure 98: Circumferential velocity in the stationary frame, with abradable grooves .............. 113 Figure 99: Swirl coefficient, with and without abradable grooves ............................................. 113 Figure 100: Meridional velocity in the stationary frame, no abradable grooves ........................ 113 Figure 101: Meridional velocity in the stationary frame, abradable grooves ............................. 114 Figure 102: Circumferential velocity in the stationary frame, surge .......................................... 114 Figure 103: Circumferential velocity in the stationary frame, midpoint .................................... 115 Figure 104: Circumferential velocity in the stationary frame, choke ......................................... 115 Figure 105: Meridional velocity, surge ....................................................................................... 115 Figure 106: Meridional velocity, midpoint ................................................................................. 115 Figure 107: Meridional velocity, choke ...................................................................................... 116 Figure 108: Velocity vectors in stationary frame, surge, 4th cavity CCW flow ......................... 116 Figure 109: Velocity vectors in stationary frame, surge, 9th cavity CCW flow ........................ 117 Figure 110: Velocity vectors in stationary frame, midpoint, 4th cavity CCW flow ................... 117 Figure 111: Velocity vectors in stationary frame, midpoint, 9th cavity CCW flow ................... 117 Figure 112:Velocity vectors in stationary frame, choke, 4th cavity CCW flow......................... 118 Figure 113: Velocity vectors in stationary frame, choke, 9th cavity CCW flow........................ 118 Figure 114: Swirl coefficient, choke to surge ............................................................................. 119 Figure 115: Direct stiffness coefficient....................................................................................... 122 Figure 116: Cross coupled stiffness coefficient .......................................................................... 122 Figure 117: Direct Damping ....................................................................................................... 123 Figure 118: Effective damping ................................................................................................... 123 Figure 119: Cross-coupled damping ........................................................................................... 124 Figure 120: Direct stiffness from choke to surge........................................................................ 125 xii Figure 121: Cross coupled stiffness from choke to surge ........................................................... 126 Figure 122: Direct damping from choke to surge ....................................................................... 126 Figure 123: Effective damping from choke to surge .................................................................. 127 Figure 124: Cross coupled damping from choke to surge .......................................................... 127 Figure 125: Circumferential velocity in the stationary frame, 20% axial groove width ........... 129 Figure 126: Circumferential velocity in the stationary frame, 40% axial groove width ............ 129 Figure 127: 3D streamlines in tooth cavity, 20% in. axial width ............................................... 129 Figure 128: 3D streamlines in tooth cavity, 40%. axial width.................................................... 130 Figure 129: Swirl coefficient, axial groove width comparison................................................... 130 Figure 130: Labyrinth tooth recirculation zone with reduced axial width, 20% in. ................... 131 Figure 131: Labyrinth tooth recirculation zone with increased axial width, 40%. ..................... 131 Figure 132: Direct stiffness, axial groove width comparison ..................................................... 132 Figure 133: Cross coupled stiffness, axial groove width comparison ........................................ 133 Figure 134: Direct damping, axial groove width comparison .................................................... 133 Figure 135: Effective damping, axial groove width comparison ................................................ 134 Figure 136: Cross coupled damping, axial groove width comparison........................................ 134 Figure 137: Circumferential velocity in the stationary frame, 20% in. radial depth .................. 135 Figure 138: Circumferential velocity in the stationary frame, 40%. radial depth ...................... 135 Figure 139: Swirl coefficient, radial groove width comparison ................................................. 136 Figure 140: Labyrinth tooth recirculation zone with reduced radial depth, 20% in. .................. 136 Figure 141: Labyrinth tooth recirculation zone with increased radial depth, 40%. .................... 137 Figure 142: Direct stiffness, radial depth comparison ................................................................ 138 Figure 143: Cross coupled stiffness, radial depth comparison ................................................... 139 Figure 144: Direct damping, radial depth comparison ............................................................... 139 Figure 145: Effective damping, radial depth comparison ........................................................... 140 Figure 146: Cross couped damping, radial depth comparison .................................................... 140 xiii KEY TO SYMBOLS A Area a major axis of rotor elliptic whirling orbit b minor axis of rotor elliptic whirling orbit C Direct damping 𝐶𝑒𝑓𝑓 Effective damping 𝐶𝑟 Radial seal clearance 𝐶𝜇 k-𝜀 turbulence model constant c Cross coupled damping 𝐷 Diameter 𝐷𝑖𝑗 Rotor displacement matrix 𝑒̃ Specific shaft work F Force 𝑓𝑜 Base frequency g Gravity constant 𝐻𝑖𝑗 Force impedance matrix 𝐻𝐼𝑠𝑒𝑛 Isentropic head K Direct stiffness k Cross coupled stiffness L Seal length M Direct inertia 𝑀̇ Mass flow rate xiv m Cross coupled inertia P Static pressure 𝑄̇ Volume flow rate q Swirl R Universal gas constant r radius T Temperature t Time 𝑢+ Dimensionless near wall velocity 𝑢𝑡 Friction velocity V Absolute velocity WFR Whirl frequency ratio 𝑦+ Dimensionless sublayer-scaled distance from the wall Z Compressibility factor Greek Letters ε Eccentricity; turbulence eddy dissipation per unit mass 𝜅 Specific heat ratio 𝜌 Density 𝜏 Time steps Ω Whirling Speed (rad/s); Excitation frequency (Hz) 𝜔 Rotor rotating speed Subscript 1 Impeller inlet xv 2 Impeller outlet n Normal 𝜃, 𝑡 Tangential x, y, z Cartesian coordinates xx,yy Direct xy,yx Cross-coupling xvi Chapter 1: Background and Motivation 1.1. Centrifugal Compressors The development of turbomachinery has a long history, dating as far back as the 1st century. [1] notes Heron in Alexandria developing a primitive steam turbine in 60 AD. Over the next millennia, turbomachinery mainly found its application in water wheels used for milling flour, and simple windmills [2]. Almost 1700 years after Heron’s steam turbine, Leonard and his son Albert Euler developed their own experiments which led to “Euler’s equation” which helped pave the way for the development of modern turbomachinery as we see it today [2]. The development of centrifugal compressors is most evidently seen in the last 80 years, paralleling closely with the advancing development of the gas turbine. An efficient working centrifugal compressor was paramount to a successful working gas turbine and the growing power needs of the 20th century. However, widespread industrial applications of centrifugal compressors began earlier around the 19th century. At this time, common uses of centrifugal compressor were used for ventilation purposes in the mining and pneumatic conveyor belts in cotton plants [2]. One of the first pioneers of the centrifugal compressor is often attributed to Auguste Rateau in the late 19th century [3]. Before then, one would fine the centrifugal impeller used in a centrifugal pump as early as 1818 when the Massachusetts Pump factory first began commercially producing them in the United States [4]. In its simplest definition, a centrifugal compressor (or pump) is a spindle wheel that adds energy to a fluid continuously, by the dynamic action of the impeller blades. In this process, the impeller imparts rotational momentum to the gas which is converted to pressure energy. The compressible working fluid density distinguishes centrifugal compressors from centrifugal pumps. This means that the fluid density in centrifugal compressors is a function of both temperature and 1 pressure. In contrast to turbines, compressors produce an adverse pressure gradient and, until recently, centrifugal pumps and compressors had low efficiencies, ranging between 45 and 60 percent. Today, centrifugal compressors have been able to reach efficiencies as high as 90 percent [4], and the demand for higher efficiencies continues to grow. 1.2. Industrial Compressors Since the 1950’s, centrifugal compressors have assumed a unique role in the natural gas and oil processing operations [5]. Compared to axial compressors, centrifugal compressors are able to achieve higher pressure ratios per stage reducing the number of stages and thereby excel when space is limited. Additionally, centrifugal compressors generally require minimal maintenance and can operate over a wide operating range with long lifecycles [5]. For the oil and gas markets, industrial compressors can be sub-divided into two categories based on their application: process and pipeline gas compressors. Process gas compressors are multi-staged compressors used for upstream applications. This includes gas gathering, gas lift, and reinjection services [5]. Pipeline centrifugal compressors are used in midstream applications such as gas transmission application. These compressors boost the pressure in pipelines, transporting process gas through miles of pipe. They are designed to accommodate high flow rates and wide operating ranges and consist of low stage count when compared to process compressors. 1.3. Introduction to Compressor Seals Without the addition of work, the fluid will flow between areas of high pressure to areas of low pressure. In centrifugal compressors, this occurs between the rotating and stationary components of the machine. Unintended admission of this fluid flow is known as leakage. Therefore, sealing systems are necessarily implemented and are critical components of centrifugal gas compressors functioning as a sealing apparatus between various sections within the system. They can dictate the efficiency, stability, and operability of the machine. Over the years, the 2 compressor industry has developed a variety of seal types, each serving a different purpose. There are principally two categories of seals: (1) shaft-end, (2) internal seals. Shaft-end seals can be further categorized as the following: (a) oil film (wet), and (b) dry gas seals [6]. This work will focus on the aerodynamic and rotordynamic performance characteristics of internal labyrinth seals. However, before moving forward, it is useful to distinguish each seal classification and their purpose, which is provided in this section. 1.3.1. Shaft-End Seals Shaft-end seals are designed to seal the rotating shaft assembly in the stationary casing. In the process industry, many applications use contacting seals or positive sealing. In this type of seal, the process gas is completely controlled, meaning that no loss of process gas can be realized. These types of seals are used when environmental regulations and safety demand that the leakage into the environment be minuscule. Thus, shaft-end seals prevent process gases from leaking through compressor components into the atmosphere. In a similar fashion, contacting seals are used to prevent atmospheric gases from mixing with and contaminating the process gas. Moreover, the process gas will usually be kept away from the compressor bearings [5]. Conventional contact sealing systems commonly used wet seals. This liquid barrier in wet seals usually comes in the form of lubricating oil [6]. The oil is circulated at high pressure around the rotating shaft, whereby the leakage gas is absorbed or blocked during the oil recirculation process. Very little gas escapes this barrier of oil and much more of it is absorbed thus mitigating leakage. In contrast to its traditional counterpart, dry gas seals have become an integral part in modern day centrifugal compressors. Favor for dry gas seals over liquid film seals is largely due to the simplicity of the system, compared to the liquid film oil system. By eliminating the complex oil lubrication system, dry gas seal applications have seen promising improvements. Operation of dry gas seals is almost 3 identical to a mechanical seal; however, it consists of a mating ring (rotating) and primary ring (stationary) controlled by hydrostatic and hydrodynamic forces which dictate the mating gap between each of the surfaces [7]. In some cases, non-contacting seals such as labyrinth seals may be used to mitigate leakage into atmosphere. To utilize non-contacting seals, a small amount of controlled leakage must be acceptable. Therefore, these seals are mainly suitable for air compression applications [4]. 1.3.2. Internal Labyrinth Seals In addition to shaft-end seals, labyrinth seals are commonly used in between compressor stages and along the balance piston in high pressure turbomachinery applications. Labyrinth seals operate between rotating and stationary parts between stages and are used to control leakage from high-pressure to low-pressure regions. The basic geometry includes a finite row of teeth located on either the rotor, stator, or both, with a small cavity in between each tooth. Typical labyrinth seal configurations are depicted in Figure 1. Figure 1: Labyrinth seal configurations a) TOS b) TOR c) Interlocking 4 1.4. Gas Labyrinth Seal Applications Common areas of utilizing of labyrinth seals in centrifugal compressors are inter-stage shaft, compressor eye, and balance piston seals [8]. These locations are illustrated in Figure 2 [9]. The function of the seal varies depending on location. Interstage shaft seals are located on the hub of the impeller and are used to inhibit flow from traveling back into the previous stage. This leakage flow is caused by the small pressure difference between consecutive stages. A second form of leakage flow occurs across the front of shrouded impellers which is produced by the pressure rise at the impeller exit causing the fluid to leak into the shroud cavity, where it recirculates back toward the impeller inlet. Compressor eye labyrinth seals are used to impede this leakage across the front of the impeller. Lastly, another form of leakage, typically of higher amount, occurs across the balance piston located at the rear of in-line compressors after the last impeller stage. The balance piston is needed to compensate for the aerodynamic thrust load enacted on the rotor due Figure 2: Typical labyrinth seal leakage for a centrifugal compressor [9] 5 to the pressure rise through the compressor [8]. Balance piston labyrinth seals are generally longer seals located on the balance piston. Moreover, the back face of the balance piston seal is connected to the compressor suction nozzle by a recirculation valve. Consequently, a percentage of the main flow at discharge leaks across the balance piston and is recirculated toward compressor suction. Therefore, balance piston labyrinth seals are critical components in reducing the flow recirculation back towards the compressor inlet. When fulfilling the objective of reducing secondary leakage flow, there is another important operating criterion associated with these seals, that must be satisfied. It is imperative that these components operate with a stable rotordynamic operating performance. The tortuous flow path inside labyrinth seals has been known to create driving forces inside the seal that can have destabilizing effect on the machine. When these destabilizing forces reach a certain magnitude, unstable vibrations can occur. Furthermore, the reaction forces associated with the rotor displacements are characterized by stiffness and damping coefficients – “rotordynamic coefficients”. Rotordynamic coefficients are used to quantify the stability of the compressor. For this reason, extensive efforts over the past decades have focused on accurately predicting rotordynamic coefficients. 1.5.Objectives The significance of studying rotordynamic behavior in gas compressor seals is a consequence of the ever-increasing demands of improving compressor performance. Designers look at every component inside the compressor in order to realize this objective. Labyrinth seals provide a way of controlling excess leakage by evolving seal designs and tightening up clearances. A reduction in seal leakage will always be accompanied by an increase in turbomachinery efficiency. However, a long history of experimental data has shown that labyrinth seals are prone 6 to induce rotordynamic instabilities at certain conditions, which can be detrimental to the machine. Therefore, more accurate predictions of labyrinth seal rotordynamic terms are desired as performance characteristics are marginalized. Simple bulk flow codes have served this purpose in compressor design over the decades with fast computational times and reasonable predictions. However, the development of CFD and advancing computational resources has largely overtaken bulk flow models in the industrial setting by providing a means of accurately predicting rotordynamic coefficients without making the same simplifying assumptions as bulk flow methods do. CFD is particularly advantageous when modelling complex geometries. The scope of this research will aim to improve the methodologies used to detail the aerodynamic fluid behavior and to predict rotordynamic coefficients found in centrifugal compressor labyrinth seals using commercial CFD software. At present, the computational resources and time required to perform unsteady analysis on an entire compressor are just becoming realistic. A single flow condition uses to require a month or more of computational running time. However, with faster CPU speeds, parallel processing, and enhanced CFD algorithms and post-processing capabilities, it is now possible to test and compare transient CFD methods for complex domains. Therefore, CFD analysis will be run using both steady and unsteady (transient) CFD methods. Moreover, the various unsteady CFD techniques will be explored, including circular and elliptical periodic orbit models in the time and frequency domains. It is predicted that the transient CFD results will provide the most accurate way of predicting rotordynamic terms. However, with its added complexity, the solution is likely to come at the cost of considerable increases in computational time. Therefore, a systematic comparison between steady-state and transient analyses will be conducted in order to determine the respective spheres of influence for each method. Should a considerable discrepancy be found between the two, it is 7 also of value to demonstrate whether or not there is a consistent discrepancy between transient and steady-state results that can be numerical reconciled. In other words, if there is a considerable difference between the two methods, is that difference consistently the same for various cases. Other considerations for this research will be examining the influence of different flow conditions on the rotordynamic characteristics of the seal. Past studies have only considered a single operating flow point, usually at the design point, while changing the geometric parameters. However, as mentioned above, process gas compressors are often needed to run at off-design points. This study will be composed of different flow conditions ranging from choke to surge in order to see how the stability characteristics behave as they go from each end of the compressor map. Additionally, many industrial gas labyrinth seals use a smooth abradable stator surface material, which allows the labyrinth teeth to rub into the abradable material. The fluid behavior and rotordynamic response created by the wear grooves in this process will be modeled and investigated in this work. 8 Chapter 2: Evolution of Theory and Flow Through Labyrinth Seals 2.1. Flow Structure and Features in the Labyrinth Seals Labyrinth seals are crucial components of centrifugal gas compressors used to control secondary leakage and recirculation loss, especially when performing at high pressure operating conditions. The flow mechanism of the labyrinth seal can be described as follows: The fluid that proceeds from the main flow path to the secondary flow path is driven by a pressure gradient which drives the fluid through the seal clearances located between the rotor and stator. The fluid continues through the tooth clearances with high velocity in a jet stream flow. Then, the streamline path expands into the cavity where its circumferential velocity increases, creating a vortex. As the fluid travels through the seal, the pressure head is converted into kinetic energy, which is dissipated by the turbulence-viscosity interaction in each passing cavity. This results in a significant pressure drop at the seal outlet and reduced leakage in the secondary flow path. Figure 3 provides an example of the flow path through a TOS labyrinth seal. Figure 3: Example of the flow path through labyrinth seal Excessive leakage (secondary flow) can represent a substantial loss in efficiency and an increase in power consumption. This is because work is being done on the recirculated fluid with no additional contribution to the compressor pressure rise. In other words, recirculation represents a parasitic loss of work that does not contribute to the specific shaft work. The specific shaft work is defined by the Euler equation: 9 𝑒̃ = 𝜔[(𝑟𝑉𝜃 )2 − (𝑟𝑉𝜃 )1 ] (1) Aungier [9] demonstrates an example of this when comparing a mean streamline performance analysis with experimental data. Moreover, in the case of the impeller eye seal, the fluid exits the impeller at a higher temperature, 𝑇2 , which then recirculates across the shroud cavity, back toward the impeller inlet. An elevated suction temperature, 𝑇1 , will affect the pressure ratio (PR) by reducing it. This is demonstrated in the relationship 𝑇1 and the PR have in the definition of isentropic head, found in Equation 2. If 𝑇1 increases, the compressibility factor, 𝑍, will increase. This is slightly counterbalanced by a decrease in the specific heat ratio, 𝑘, however, the PR will also change. Therefore, if the same 𝐻𝑖𝑠𝑒𝑛 is maintained, the pressure ratio will decrease. An example of this is found in the CFD work of Qiao et al. [10] which investigated the effects of leakage on a multi-stage centrifugal compressor. Qiao et al. [10] found that when modelling the secondary leakage flow, the pressure ratio was reduced when compared to a model without leakage. In a similar relationship, if the same PR were to be preserved, a higher 𝐻𝑖𝑠𝑒𝑛 is needed, which requires a faster running speed and, thereby, more power. 𝜅−1 𝜅 𝑃2 (2) 𝜅 𝐻𝐼𝑠𝑒𝑛 = 𝑅𝑍𝑇1 ( ) [(𝑃 ) − 1] 𝜅−1 1 Additionally, multi-stage compressors experience leakage across the hub of the rotor through the impeller shaft labyrinth seals. At this location, a portion of the flow that is traveling through the return guide vane to the next compressor stage, circulates back through the shaft labyrinth seal toward the preceding stage impeller exit. This recirculation increases the turbulence due to mixing. The flow mixing loss contributes to a reduction in compressor efficiency [11]. A third form a of leakage is experienced in inline multi-stage compressors after the last stage, through the balance piston labyrinth seal. Unlike the previous seals mention, the balance 10 piston labyrinth seal experiences pressure differentials of the entire compressor. The pressure upstream of the balance piston seal is approximately the pressure leaving the last stage impeller. The flow downstream the balance piston seal goes back to compressor suction, creating a known pressure ratio. Moreover, the fluid that is recirculated back toward suction is worked on again through each stage of the compressor, consuming more power. For this reason, the leakage through the balance piston labyrinth seal has the greatest influence on compressor efficiency. However, in addition to this, the added leakage mass flow combined with a higher inlet temperature contributes to a higher inlet volumetric flow rate. This relation is described in Equation 3. The high-volume flow increases the surge margin by pushing the operating point further out on the compressor map, diminishing the head [8]. 𝑀̇𝑅𝑇 (3) 𝑄̇ = 𝑃 Up to this point, it has been demonstrated that annular gas labyrinth seals play a vital role in maintaining desired operating efficiencies by control leakage losses. One of the best ways to reduce the leakage in these seals is to minimize the clearances. However, this is typically easier said than done. At high operating speeds, the rotor will experience centrifugal growth, and therefore, labyrinth seal teeth with tight clearances run the risk of breaching the clearance gap, causing contact between the rotor and stator (rubs), thereby inducing unstable vibrations. Improvements in this area have come in the development of abradable seals. This solution implements a softer abradable surface on the stator that allows for the TOR labyrinth teeth to cut into the surface if contact should occur [12]. This allows for tighter tooth clearances without experiencing the detrimental effects of rotational rubs. Other geometrical parameters effect leakage as well, such as tooth shape, pitch/depth, and number of teeth. Operating conditions (i.e., temperature, shaft speed, supply pressure and discharge pressure) and gas type are also important 11 contributing factors to seal leakage [13]. Additionally, the configuration of the flow path of the recirculating fluid exiting the impeller eye seal must be designed to accommodate a smooth transition from the secondary flow path back into the main flow path [3]. Figure 4a. represents a poor design as the recirculating flow path creates a vortex upstream of the impeller inlet, diminishing the compressor efficiency [3]. Figure 4: shroud face configurations of impeller eye seal a) open face, b) straight face c) curved face. [3] 2.2. Rotordynamics Rotordynamics is a branch of engineering that deals with applying the general laws of dynamics to rotating systems and their surrounding structures. In the realm of turbomachinery, this involves the study of lateral and torsional vibrations of rotating shafts with the objective of predicting and controlling those vibrations within the system [13]. During the operation of a centrifugal compressor, these vibrations can be caused by two kinds of mechanisms: mechanical excitations and aerodynamic excitations. Mechanical excitations include rotor imbalances or rotor rubs in areas of tight clearances. Rotor imbalance occurs when there is an unequal distribution of mass about the rotor centerline. Rubs occurs when part of the rotor contacts the stator. Aerodynamic excitations are caused by the fluid-structure interactions within the various components of the 12 compressor. Aerodynamic excitation can occur when the impeller or diffuser experiences rotating stall or if the compressor goes into surge. Moreover, the occurrence of aerodynamic excitation can also be caused by the leakage flow path through the labyrinth seals. It is this form of aerodynamic excitation that will be the focus of this work. The classical model used to explain the basic concepts of rotordynamic systems is referred to as the Jeffcott rotor. The Jeffcott rotor illustrates a single mass rotor (disk) centrally located on a massless shaft supported by two bearings (Figure 5). In this idealized case, a flexible rotor with rigid bearings is considered. Though simple, this system gives valuable insight to the mechanics of rotor whirl due to mass imbalance and provides the framework for understanding more complex systems. Mass imbalance causes the center of gravity of the rotor to become eccentric to the undisturbed axis of rotation. A rotor is said to whirl when the center of gravity of any cross section traces out an orbit in time, instead of remaining at a fixed point [14]. An example of this whirling motion is depicted in Figure 6. The whirling rotor problem that the Jeffcott rotor addresses is if the whirl motion is in the same direction as the rotational velocity, it is classified as forward whirl. Conversely, if it is in the opposite direction of the rotational velocity, it is considered backward whirl. In Jeffcott’s simple model, he explained how the rotor whirl amplitude is maximized at the critical speed, that is, the speed which coincides with the systems natural frequency, and begins to diminish beyond that [15]. 13 Figure 5: The Jeffcott rotor [18] Figure 6: Whirling Rotor Model [16] 14 2.2.1. Rotordynamic Coefficients When the fluid travels through the gas labyrinth seals, the circumferential flow path of the fluid imposes net pressure and shear forces that act upon the rotor. The reaction forces induced by the rotor displacements are quantified by direct and cross-coupled stiffness and damping coefficients, otherwise known as rotordynamic coefficients. Figure 7 displays a rotordynamic model for a labyrinth seal with a centered rotor, and the relationship of the reaction forces acting on a whirling rotor are shown in Figure 8. Corresponding to these two figures, a brief description of each rotordynamic coefficient is provided in this below. Moreover, for simplicity, the following description corresponds to a circular whirl orbit. Figure 7: Rotordynamic model of gas labyrinth seal with a centered rotor 15 Figure 8: Direct and cross-coupled seal reaction forces acting on a whirling rotor The direct stiffness coefficient as shown in Figure 7 is denoted as 𝐾𝑥𝑥 = 𝐾𝑦𝑦 . In General, the direct stiffness values for labyrinth seals are low and do not have a major influence on stability. But in principle, when considering a whirling rotor, a positive direct stiffness value acts to center the rotor and resist radial motion. Negative direct stiffness values, which are common for labyrinth seals, produce a radially force directed outward which acts to pull the rotor radially outward toward the seal [16]. Direct stiffness values have an important effect on rotordynamic performance if a sufficiently large negative value is realized. Large negative stiffness values can indirectly affect rotor stability by influencing the vibrational frequencies that affect the rotor's critical speed [16, 17]. Additionally, large negative values have been noticed in back-to-back compressor configurations [18]. Moreover, Direct stiffness depends strongly on the number of cavities of the labyrinth seals [18]. Therefore, the behavior of this coefficient may be more important in longer labyrinth seals such as the balance piston seal. Similar to the direct stiffness term, the cross-coupled damping term, 𝑐𝑥𝑦 = −𝑐𝑦𝑥 does not have a direct effect on rotordynamic stability and has little impact on gas labyrinth seals. For this 16 reason, it is often neglected when considering rotordynamic stability. But in principle, this term produces radial force colinear with the deflection vector [16]. A positive 𝑐𝑥𝑦 value corresponds to stiffening. According to [19], the cross-coupled stiffness term is denoted as 𝑘𝑥𝑦 = −𝑘𝑦𝑥 . Due to its destabilizing influence, tends to have the greatest influence on the rotordynamic performance. When 𝑘𝑥𝑦 is positive, a reaction force in the same direction of positive whirl is realized, which promotes instability. Contrarily, a −𝑘𝑥𝑦 term tends to promote rotor stability. Cross coupling forces become more significant in high pressure and rotating speed applications, and when clearances between the rotating and stationary components are small [20]. Additionally, the circumferential flow – "swirl," is often associated with the cross-couple stiffness term. Therefore, operations with higher swirl conditions tend to produce higher cross-coupled stiffness values and promote instability. For this reason, it has become common practice in industry to implement anti- swirl technology to reduce the inlet swirl entering the labyrinth seal. One common solution to reduce swirl is to implement Anti-Swirl Vanes (ASV) or swirl breaks. ASV’s have been utilized in many seal applications and are implemented at the seal entrance [20]. Another method for reducing swirl is by using a shunt injection. Shunt injections are often found in longer labyrinth seals such as balance pistons. A shunt injection directs part of the discharge pressure back towards the rotor near the high-pressure end of the labyrinth seal. Consequently, circumferential component of the flow is converted into a more axial one. To illustrate this, Figure 9 and Figure 10 show a balance piston seal configuration without and with a shunt injection, respectively. 17 Figure 9: Conventional balance piston labyrinth seal flow [21]. Figure 10: Balance piston labyrinth seal with shunt injection [21]. The destabilizing forces are counterbalanced by the compressors damping forces. Much of the systems damping will come from the bearings. However, the seals can also provide enough damping that makes them vital components in determining stability margins [18]. The direct damping term, denoted as 𝐶𝑥𝑥 = 𝐶𝑦𝑦 , characterizes this stabilizing component of gas labyrinth seals. When positive, this term acts as drag force in the opposite direction of rotor whirl, effectively 18 counteracting the destabilizing 𝑘𝑥𝑦 term. Pelleti and Childs [18] observed the direct damping coefficient increasing with decreasing pressure ratio. This also corresponds to an increase in direct damping as the average density decreases. However, observation in Arthurs experiment [22] did not show a significant impact on of the pressure ratio on the direct damping for the TOR labyrinth seal. Moreover [22, 23, 24] observe an increase in direct damping as the inlet swirl increases in each of their experiments conducted at Texas A&M. Additionally, stability of the gas labyrinth seal is often characterized by using nondimensional parameters such as the effective damping (Ceff) and whirl frequency ratio (WFR). These parameters are defined in Equations 4 and 5. Both of these parameters relate the stabilizing direct damping term with the destabilizing cross coupled stiffness term. For Ceff, it is desirable to have a large positive value. A negative Ceff value shows that the destabilizing force exceeds the damping forces within the seal. For WFR, stability tends to decrease at larger ratios. 𝐶𝑒𝑓𝑓 = 𝐶 − Ω 𝑘 (4) 𝑊𝐹𝑅 = 𝐶Ω 𝑘 (5) 2.3. Literature Review Over the years, rotordynamic prediction methods have made great strides through experimental data, bulk flow code formulation, and in more recent years, computational fluid dynamic predictions. Experimental setups became the dawn of turbomachinery seal studies dating as early as 1965. Alford [25] used experimental methods to study the effects of self-excited rotor whirl associated with labyrinth seals inside axial compressors and turbines of aircrafts. He investigated the circumferential variation of static pressure acting on the cylindrical surface of the rotor inside the labyrinth seal as a primary origin for disturbing forces within the rotor. He also suggested the significance of the eccentricity of the rotor as a causation of perturbations on the 19 rotor. His set up involved a single cavity labyrinth seal with different flow areas at the inlet and discharge tooth. He posited that the clearance gaps of labyrinth seals at the inlet and outlet contribute to the stability of the rotor. Specifically, he noted that a minimum flow area at the outlet produces whirling in the same direction as rotation. Contrariwise, minimum flow area at the inlet produces stable conditions. Alford also proposed that to prevent the whirling problem, a more rigid rotor system could be implemented. His work would pave the way for future rotordynamic interests in labyrinth seals. Adding to the experimental realm, Bencker and Wachter [26] initiated the analysis of rotordynamic coefficients of labyrinth seals associated with different operating conditions. Counter to previous work, Bencker and Wachter studied multi-cavity labyrinth seals, verses single cavity, with a static eccentricity of the rotor. The authors found that lateral forces inside the labyrinth seal were caused by circumferential velocity of the fluid flow inside the cavities. They proposed that the reason for such velocity components had potentially two causes: One being that the drag induced by the rotating shaft produces a circumferential velocity in the cavity. The other being the boundary condition at the inlet where pre-swirl was present. The primary purpose of these experiments was to calculate spring coefficients that would be able to predict tangential forces inside the labyrinth seal. R.G. Kirk [27] of Virginia Polytech did some work investigating the current experimental theories, as of 1988, of calculating rotordynamic stability parameters and discusses the practicality of applying said theories to centrifugal compressor designs [27]. His approach was to compare various experimental design calculations to test data from well-established centrifugal gas compressors in the field. Kirk specifically examined the theories of Alford and Bencker and Wachter, mentioned above, as well as others. Alfords procedure for calculating turbomachinery 20 excitation continued to be the widely accepted approach for both axial and centrifugal turbines and compressors. According to Kirk, its use was proven to be a reliable indicator of rotor excitation, however it relied heavily upon a correct correlation factor. Alfords model also neglected circumferential flow inside the labyrinth seal. Becker’s approach, according to Kirk, proved to be useful in theory, however it became apparent that seal damping coefficients also needed to be considered along with the spring constants. Therefore, certain adjustments needed to be made to address the seal damping, to provide meaningful results. Pelletti [28] presented comprehensive experimental results for rotordynamic coefficients for a short TOS and TOR labyrinth seal. The experimental results explored the effects of inlet pressure, pressure ratio, rotor speed, seal clearance, and inlet pre-swirl. According to Pelletti’s study, increasing the rotor speed and decreasing the pressure ratio had a stabilizing effect on both seal configurations. The experimental trends were compared to the theoretical predictions of a two- control volume compressible flow model presented by Sharrer [29]. The theoretical model correctly predicted the trends associated with pressure ratio, rotor speed, and pre-swirl. However, the theory failed to correctly predict the effects of seal clearance on rotordynamic coefficients. Pelletti also suggested that since stability is enhanced through reduction in tangential velocity, incorporating an effective swirl brake upstream the seal may be beneficial in increasing stability. Kwanka [17], reports experimental data for a 5-tooth stepped labyrinth seal configuration. In this experiment, he compares the efficacy of controlling leakage using a stepped tooth configuration compared with a straight through configuration. As predicted, the stepped seal configuration leaks less. However, it is noticeable that the rotordynamic properties of the of the stepped seal are less favorable. 21 Expanding on experimental methods, the 1980s brought forth analytical methods with the use of bulk flow models to predict rotor instabilities. The development of bulk flow models proved to be effective in obtaining fast calculations. Iwatsubo [30] is noted for providing one of the first one-control-volume models using bulk flow calculations to predict rotordynamic coefficients. Iwatsubo used a finite difference approach to solve for the circumferential momentum and continuity equations which he used to define a bulk flow circumferential velocity inside the labyrinth seal. The velocity and displacement of the rotor characterized the force inside the labyrinth seal. The spring and damping coefficients related to the force and effecting rotor stability could then be described using energy conservation. He concluded that the force imposed by the labyrinth seal always makes the rotor system unstable. Wyssmann et al. [31], would later establish a multi-control volume model to calculate the labyrinth coefficients in straight labyrinth seals. His model calculated the circumferential pressure distribution as a function of lateral rotor displacement and displacement time rate. The flow configuration in his model is drastically simplified so that the time-averaged Navier Stokes equations using k-e turbulence model could be solved using a computer program. Wyssmann confirmed the importance of reducing the inlet swirl to eliminate cross-coupled stiffness, and found that when reducing the inlet swirl, the damping capacity of the seal is retained. Childs and Scharrer [32] sought to improve the single volume control model of Iwatsubo and expand upon the two-control-volume approach of Wyssmann et al. [31] by incorporating the recirculation velocity in the two-control-volume model. They compared their theory to experimental data consisting of both teeth and stator and teeth on rotor labyrinth seal configurations. Their model gave reasonable predictions for cross-coupled stiffness for both TOR and TOS see through labyrinth seals, however, failed to predict the trend of decreasing direct 22 damping coefficient with increasing clearances for the TOS configuration which was represented in experimental data. B.P Williams and R.D Flack [33] would later modify the Iwatsubo bulk flow method which they used to establish their own approach in solving rotor dynamic coefficients. Specifically, they iteratively solved for the mass flow based on pressure drops across the teeth of the labyrinth seal, and then applied that mass flow rate into the governing Navier Stokes equations. However, small variations in mass flow rate values did not affect the dynamic coefficients all that much. When comparing the bulk flow code to experimental data, the trends established proved to be useful for design criteria. The accuracy however was limited, especially with direct stiffness and direct damping. An attempt to achieve a better correlation with experimental data using a correlation parameter was made without success. Experimental and analytical analysis of labyrinth seal induced instabilities proved to be useful for many years. However, over the years, as design parameters of centrifugal compressors have been increasingly becoming more demanding, an even greater knowledge of critical design parameters associated with rotordynamic instabilities has become more necessary. Over the recent decades, computational fluid dynamics has evolved into an important method in achieving highly detailed analysis in centrifugal compressor design and analysis. Computational fluid dynamics became relevant in studying gas seals as early as the 1990’s. Earlier CFD models manifested themselves in self-developed codes which solved the 3D Reynolds averaged Navier Stokes equations using finite differencing scheme. Since then, commercial CFD codes such as Ansys have become an integral part of industrial analysis. One of the early pioneers in this type of CFD analysis was Rhode et al. [34], who developed a 3-D finite difference approach to analyze an eccentric whirling labyrinth seal. Due to the 23 computational demands of his model, only a single cavity the multi-cavity seal was run with each computation. As a result, he could calculate rotordynamic forces on the labyrinth seal and compare his numerical results with measured data. He found that the tangential force agreed much closer with measurements than the radial force. He was also able to show that these forces were increased as the inlet swirl increased. Moore [35] used a CFD code, SCISEAL, adapted from [36], to predict the flow and rotordynamic coefficients of an 8 tooth TOS labyrinth seal presented in [28], which he benchmarked with Pelletti’s experimental data, and further compared the CFD results to bulk flow predictions. He concluded that CFD was advantageous in achieving a better prediction of rotordynamic coefficients over bulk flow codes when compared to experimental data. Despite this, Moore noted that due to computational demands of CFD, bulk flow codes remain a useful tool in achieving fast analysis. Hirano et al. [37] analyzed a typical labyrinth seal of a steam turbine engine and a centrifugal compressor eye seal using a CFD code, CFX-TASCFlow. A 3D model with an eccentric rotor was used to solve for the leakage and rotordynamic forces using the CFD code. Each case was solved using an inlet pre-swirl ratio of 0. The results were then compared with existing bulk flow code DYNLAB. The comparison showed that the bulk flow models gave pessimistic predictions of the destabilizing forces compared to CFD. However, bulk flow codes were said to be advantageous at the current time for their fast-computational times. Kirk and Guo [38] discussed in another paper studying the same seal that it is possible to calibrate the friction factors in bulk flow models to get a closer match with CFD. Thompson [39] continued CFD efforts in analyzing labyrinth seal rotordynamic characteristics using commercial CFD Ansys-CFX. She analyzed a compressor eye straight 24 labyrinth seal and compared her results with a bulk flow code, VT-FAST, and another commercial CFD code, TASCFlow. The resultant forces calculated through each method showed to be quite different. The author noted the need to perform a mesh-sensitivity study in order to validate the accuracy of the results. Wen-sheng et al. [40] investigated the leakage and stability of a straight labyrinth seal (TOS) using CFD software. The stability study was conducted by comparing the rotordynamic coefficients with varying seal clearances. [40] found that as the seal gap increased, each rotordynamic term decreased. Gao [41] did a comprehensive analysis on three types of labyrinth seals: straight compressor eye (TOS and TOR), a stepped eye seal, and a balance piston seal. He modeled each seal using ANSYS-CFX. Using the straight seal, he investigated the effects of eccentricity, clearance, tooth configuration (TOS vs. TOR), and varying pre-swirl rates. He found that increasing pre-swirl rate increased the cross-coupled and direct stiffness as well as the direct damping, while the cross-coupled damping decreases with increased pre-swirl. He also concluded that direct stiffness always increases with clearance while the other coefficients either increase or decrease depending upon pre-swirl rate. With regards to the seal configuration, he found that the TOR had a larger cross-coupled stiffness (destabilizing) term but also a larger direct damping (stabilizing) term. The stepped seal was compared with a bulk flow code and experimental data. Each method for the stepped seal showed close comparisons, particularly at lower pre-swirl rates. Lastly, his analysis on the balance piston seal confirmed the trend of increasing instability with higher preswirl rates. Other CFD studies include Subramanian, who added the impact of centrifugal growth on the rotordynamic characteristics of labyrinth seals [42]. 25 Most rotordynamic CFD analysis found in literature up until now has used steady state CFD in order to make the problem solvable in a realistic computational time. Now, many companies and universities have access to powerful computational resources that are able to run transient CFD in a reasonable turnaround time. Wu et al. [43] studied the leakage and rotordynamic coefficients of an interlocking seal and straight TOS labyrinth seal using a transient CFD approach. The results showed that the ILS seal has considerably less leakage than the TOS labyrinth seal. However, for a wide range of frequencies, the TOS labyrinth seal has a higher effective damping (stabilizing) than the ILS. The authors also compared the CFD results to the bulk flow models and found a significant difference in results and therefore proposed caution when using simplified bulk flow models. Some of the most recent works include Chochua [44]and Yan et al. [45] who use a transient CFD model to predict rotordynamic coefficients of a hole pattern seal. In contrast to [44] which uses a unidirectional perturbation method, [45] gives a two-directional motion of the rotor surface by considering a concentric whirling model. Li et al. [46] used a multi-frequency elliptical whirl model to predict rotordynamic coefficients of 3 types of annular gas seals. This work demonstrated good accuracy using Transient CFD for each labyrinth seal configuration. In an additional work, Li et al. [47] compared single- frequency and multi-frequency transient whirling rotor models to predict rotordynamic coefficients of a pocket damper seal. [47] found good agreement with experimental data and each method gave comparable results to one another. However, [47] found the multi-frequency transient whirling rotor models were advantageous in faster computational times over the single-frequency models. 26 Chapter 3: Methodology 3.1. Experimental Benchmark Although experimental setups contribute as some of the oldest methods used to predict rotordynamic behavior in gas labyrinth seals, available literature identifying rotordynamic coefficients from experiment remain relatively scarce. Kwanka [48], notes one reason is that the forces introduced by compressible fluids are generally an order of magnitude smaller than liquid labyrinth seals, making them harder to measure. Therefore, experimental setups can be both time consuming and costly. Nevertheless, they provide valuable information on rotordynamic coefficients and trends associated with gas labyrinth seals. For this reason, experimental investigations of rotordynamic characteristics of labyrinth seals readily available in the literature are used to benchmark the CFD methodologies carried out in the proceeding sections. In general, each rotordynamic identification method consists of exciting the rotor (or stator) with a certain force and measuring the relative rotor-stator displacement. A time-domain or frequency-domain approach is then carried out on the seal which is modeled as a linear system, characterized by stiffness, damping, and inertia parameters. The force-motion relationship (for small motion) using a linear system of stiffness, damping, inertia terms is modeled in Equation 6 [18]. 𝐹𝑥 𝐾𝑥𝑥 𝑘𝑥𝑦 𝑥 𝐶𝑥𝑥 𝑐𝑥𝑦 𝑥̇ 𝑀𝑥𝑥 𝑚𝑥𝑦 𝑥̈ (6) -[𝐹 ] = [ ] {𝑦 } + [ ]{ } + [ ]{ } 𝑦 𝑘𝑦𝑥 𝐾𝑦𝑦 𝑐𝑦𝑥 𝐶𝑦𝑦 𝑦̇ 𝑚𝑦𝑥 𝑀𝑦𝑦 𝑦̈ For annular gas seals, the added mass terms are negligible. Additionally, for small motion about a centered position, 𝐾𝑥𝑥 = 𝐾𝑦𝑦 = 𝐾, 𝑘𝑥𝑦 = −𝑘𝑦𝑥 = 𝑘, 𝐶𝑥𝑥 = 𝐶𝑦𝑦 = 𝐶, 𝑎𝑛𝑑 𝑐𝑥𝑦 = −𝑐𝑦𝑥 = 𝐶 and the linear system is often reduced to the following [18]: 27 𝐹𝑥 𝐾 𝑘 𝑥 𝐶 𝑐 𝑥̇ (7) [𝐹 ] = [ ] {𝑦 } + [ ]{ } 𝑦 −𝑘 𝐾 −𝑐 𝐶 𝑦̇ Moreover, the rotordynamic coefficients for gas labyrinth seals are typically frequency- independent. However, some recent experimental investigations have evidenced frequency- dependent characteristics, particularly at higher excitation frequencies [19]. Wagner et al. [49] observed that the frequency-dependent term in these cases is predominantly in the radial dynamic stiffness coefficients. Moreover, Wagner et al. [49] notes that the quadratic frequency-dependency radial term is accounted for by the direct inertia term in Equation 6. Childs [18] presents a wide collection of experimental identification of rotordynamic coefficients in gas labyrinth seals in his book “Turbomachinery Rotordynamics”. One study references an experimental investigation of a short labyrinth seal conducted by Joseph Pelleti at Texas A&M [28]. This study presents the first comprehensive study of short labyrinth seals that measures both stiffness and damping coefficients. For this reason, the fluid dynamic and rotordynamic characteristic of the experimental seal is used to validate the CFD methodologies in this work. To identify the rotordynamic coefficients in Pelleti’s experiment, the rotor is excited by a hydraulic shaker. The hydraulic shaker is used to predefine the perturbation displacement and the accompanying force response from the prescribed motion allows for the rotordynamic coefficients to be measured. This identification process is described in detail in Boettler et al. [50]. Additionally, this test set up also measures the axial pressure distribution and leakage characteristics, which are used to validate the CFD aerodynamic predictions of the seal. Pre-swirl is introduced upstream the labyrinth seal inlet via inlet swirl guide vanes. An additional benchmarking case using a hydraulic shaker is taken from the experiment set up by Arthur [22]. 28 This experiment supplements benchmarking data by offering higher pressure conditions and a longer TOR labyrinth seal compared to the Pelletti experiment. Another method of identifying rotordynamic coefficients in gas labyrinth seals is through the use of magnetic bearings. This method, inspired by Ulbrich [51] , is used in the experimental investigation of rotordynamic coefficients carried out by Kwanka [17]. The magnetic bearing is used to support the rotor and also serves as a force transducer which allows for the reaction forces in the tangential and radial direction to be measured directly. The testing in this experiment was performed under low speed and low-pressure conditions making it unrealistic for common industrial applications. Nevertheless, the test data and established rotordynamic trends make it a suitable benchmarking choice to compare with the CFD. 3.2. Computational Fluid Dynamics The primary goal of this study is to construct a well-validated CFD methodology to predict aerodynamic performance and to quantify the rotordynamic characteristics of gas labyrinth seals used in high performance centrifugal gas compressors. To fulfill this objective, the commercial CFD software ANSYS CFX will be used to run the simulations described in this work. This software was chosen because it is a readily available and commonly used in the process industry and it easily accessible in the academic setting. Additionally, it provides the choice of several turbulence models suitable for gas compressor applications. The two most frequently used turbulence models in turbomachinery applications are the 𝑘-𝜀 and 𝑘-𝜔 Shear Stress Transport (SST) model, which will be considered and compared and are discussed in further detail in section 3.4. The complex modeling of the labyrinth seal fluid domain requires the use of high- performance computers to run each simulation. This work takes advantage of a set of parallel 29 processing high performance computers provided by the Institute for Cyber-Enabled Research at Michigan State University. Until recently, transient CFD simulations performed on labyrinth seals have been hindered by unrealistic computational times. Parallel processing has matured to a point where unsteady CFD can be computed in reasonable turnaround times. This has made it possible to compare steady-state and unsteady CFD methodologies, which will be used to draw conclusions on whether there is an added benefit in using a more computationally expensive transient CFD approach. 3.3. Computational Models 3.3.1. Tooth on Rotor Labyrinth Seal A set of TOR labyrinth seal configurations are used to study the fluid dynamic and rotordynamic performance – two of which are found in the experimental literature. The first seal analyzed, referenced in the literature, is an 8-tooth on rotor straight labyrinth seal tested by Pelletti [28]. Seal 4 in [28] and the corresponding boundary conditions are used for this analysis. One can also find the results of this study in Child’s book [18]. However, when specifying the pressure boundary conditions for this seal according to boundary conditions prescribed in [28], it was found that the leakage value was substantially higher than that reported in the experiment. Moore [35] reported a similar problem when using CFD to analyze the 8 TOS (seal 1) from Pelletti's experiment. Moore suggested that the clearances reported in [28] might not have been confirmed by measurement. Following this advice, the clearance used in the presented CFD model was reduced until the predicted CFD leakage correlated with experiment. The modified clearance is 0.105 mm. This correction is validated in Figure 11, which shows that with the corrected clearance, the predicted pressure drop across each cavity matches well with the experimental pressure measurements in each cavity. 30 Figure 11: Comparison of predicted pressure distribution with experimental measurements In a more recent work, Arthur [22] carried out experiments on a TOR labyrinth seal with high operating pressures. In this work, the leakage and rotordynamic coefficients are measured at 3 rotational speeds, pressure ratios, and pre swirl ratios. Furthermore, this experiment provides investigation of a long labyrinth seal with 14 teeth and a smooth stator. Therefore, this seal is chosen to validate further CFD investigation. The dimensions and operating conditions for the experimental labyrinth seals presented by Pelletti and Arthur are summarized in Table 1. Table 1: Boundary conditions for TOR labyrinth seals 8-Tooth on rotor [28] 14-Tooth on rotor [22] Number of teeth 8 14 Seal radius 75.679 mm 57.15 Clearance 0.105 mm 0.10 mm 𝐿 𝐷 0.166 0.620 Working fluid Air (ideal gas) Air (ideal gas) Inlet pressure 18.3 bar 70 bar Inlet temperature 306 K 288 K Pressure ratio 0.66 0.5 Inlet swirl ratio 0, 0.165 0.08, 0.25, 0.48 Rotating walls 1675.52 rad/s 1068.14 rad/s 31 Lastly, a CFD analysis performed on the balance piston labyrinth seal belonging to Solar Turbines Inc. C31 compressor. This seal is characterized by 18 teeth and will be analyzed under different flow conditions ranging from choke to surge in order to see how the stability characteristics behave as they go from each end of the compressor map. The flow points on the Speedline are contributed by Solar Turbines. An example of the flow conditions along the compressor map are shown in Figure 12. Moreover, the balance piston labyrinth seal operates with a shunt injection called P2 inject. The boundary conditions will be applied in accordance with Figure 13. For simplicity, P2 inject is modeled as a circumferential slot instead of discrete holes located around the circumference of the seal. 70000 60000 Isentropic Head (FT-LBf/LBm) 50000 40000 15300 rpm 30000 cfd points 20000 10000 600 800 1000 1200 1400 1600 1800 2000 Inlet Volume Flow (CFM) Figure 12: Operating flow points along the compressor map 32 Figure 13: Boundary Conditions Used for the Balance Piston Simulations Table 2: Boundary conditions, balance piston labyrinth seal Pressure (psi) Temperature (F) Condition RPM P2 Inject Opening Outlet P2 Inject Inlet Inlet Inlet Choke 15300 3351 3217 1320 254 Midpoint 15300 3985 3786 1320 283 Surge 15300 4170 3987 1320 294 3.3.2. Tooth on Stator Labyrinth Seal In addition to the TOR configuration, a set of TOS labyrinth seals are also investigated. The first seal described here is a straight 5-tooth on stator seal. This prototype seal was first introduced by Hirano et al [37] and has since been the focus of previous CFD research found in the literature [41, 42]. The purpose of this analysis is to determine the reliability of steady-state CFD models in comparison to other CFD and bulk flow models when experimental data is not available. Next, a 5-tooth stepped labyrinth seal presented by Kwanka [17] in his experimental work, is also examined. Gao et al. [52] also performed a CFD analysis on the same stepped labyrinth seal. The CFD results published in this work are also included for comparison. The stepped-tooth design is an attempt to improve leakage control when compared to a conventional straight labyrinth seal. For this study, it is noticeable that the stepped eye seal has much larger pre- 33 swirl ratios. In most cases, well over 100%. This is because the rotating speed of the stepped eye seal is much lower compared to normal industrial conditions, making it possible to achieve much higher swirl ratios. In reality, swirl ratios typically found in high-speed industrial labyrinth seals are less than 100%. However, since experimental data was provided for this seal, these conditions were used to validate the CFD. Table 3: Boundary conditions for TOS labyrinth seals 5-Stepped Eye Seal [17] 5-Tooth on stator [37] Number of teeth 5 5 Seal radius 90 mm 137.4 mm Clearance 0.5 mm 0.292 mm Working fluid Air (ideal gas) Air (ideal gas) Inlet pressure 0.201 MPa 3.447 Mpa Inlet temperature 300 K 366.7 K Pressure ratio 0.502 0.5 Inlet swirl ratio 0,4,8,10,13,19 0, 0.468, 0.736 Rotating walls 78.54 rad/s 1162 rad/s 3.4. Turbulence Models As mentioned above, this work uses Ansys CFX to solve for the discretized Reynolds Averaged Naiver-Stokes (RANS) equations. The additional unknowns produced during the averaging procedure – i.e., Reynolds Stresses, are modeled using turbulent models available in CFX. The most common turbulence models used in turbomachinery applications are the k-ε and the Shear-Stress-Transport k-ω based model. These models are based on the eddy viscosity hypothesis which relate the Reynold stress to the mean strain rate. The k-ε model is the most widely used turbulence model in CFD, mainly due to its robust characteristics and a broad range of applicability. In this solver, the turbulent fluctuations are related to the mean velocities through a known turbulent viscosity. Furthermore, two transport equations are solved for the two turbulent quantities: the turbulent kinetic energy, k, and the turbulent eddy dissipation, ε. To improve accuracy, a high-resolution advection scheme is used, 34 which specifies a blending factor, between 0 and 1, and pushes the solver towards a second order upwind differencing scheme. Furthermore, this model uses five constants based on empirical observation to give an accurate estimation for a variation of flows [53]. Based on the eddy viscosity concept, the turbulent kinetic energy (k) and the dissipation rate (𝜀) are related to the turbulent viscosity in the following way: 𝑘2 (8) 𝑢𝜏 = 𝜌𝐶𝜇 𝜀 The 𝑘-𝜀 model is valid for modelling fully turbulent flows but performs poorly in the boundary layer region near the wall. ANSYS CFX, accommodates for this by using a wall function approach. Wall functions are empirical functions fitted to the observed universal behavior of what is known as the Law of the Wall (Figure 14). The law of the wall splits the turbulence boundary layer into three regions: the viscous sublayer, the buffer layer, and the log-law region, which relate the average velocity of a turbulent flow to the distance from that point to the wall. Furthermore, this relationship defines the near wall velocity parallel to the wall and wall distance normal to the wall in dimensionless terms of 𝑢+ and 𝑦 + , respectively. A third important term is the friction velocity, 𝑢𝜏 , which characterizes the shear stress in terms of velocity. Moreover, 𝑢𝜏 completes the 𝑢+ definition given below. 𝜏 (9) 𝑢𝜏 = √ 𝜌𝜔 𝑢 𝑢+ = 𝑢 (10) 𝜏 𝑦𝑢𝜏 𝑦+ = (11) 𝜈 Referring to Figure 14, the region of 𝑦 + < 5 is referred to as the viscous sublayer. As the name suggest, this region is dominated by the viscous forces, and the flow can be considered 35 essentially laminar [54]. Furthermore, the Reynolds Stresses are negligible, leaving a linear velocity relationship in Equation 13. 𝑢+ = 𝑦 + 13) At 30 < 𝑦 + < 300, the turbulent stress dominates, and the velocity varies with a logarithmic function along the distance 𝑦 + , distinguishing this region as the logarithmic region (log-wall). 1 (14) 𝑢+ = ln (𝑦 + ) + 𝐶 𝑘 Figure 14: Law of the Wall 𝑘 is the von Karmen constant and C is the log-wall constant based on wall roughness. Adjoining the viscous sublayer and the log-wall region is the buffer layer, 5 < 𝑦 + < 30. In this region, the molecular viscosity and turbulence are of equal importance and the velocity profile is not as clearly defined. Most CFD codes implement a proprietary blending function of the viscous sublayer and logarithmic region to capture the effects of the buffer layer. 36 Using the definitions above, the wall functions approach attempts to bridge the viscosity effected sublayer and the fully-turbulent region using the empirical formulations. Without wall functions, one would have to resolve the boundary layer by integrating the turbulence to the wall. This would require a very fine mesh in order to ensure that the first cell center is placed in the viscous sublayer (𝑦 + ≈ 1). However, this is difficult to achieve and requires additional viscous damping functions to avoid numerical instabilities and singularities near the wall [55]. To avoid this, ANSYS CFX implements scalable wall functions to prevent problems of successive refinements in standard wall function meshes. Scalable wall functions ensure that the wall distance applied is such that y+ ≥ 11.06 For very fine grids, the wall shear stress is limited such that y+ = 11.06. For grids courser than y+ = 11.06, the results are obtained using standard wall functions [53]. The advantage of this model is the ability to use a courser mesh and achieve faster computational results. CFX also has the option of implementing a k-ω Shear-Stress-Transport model. In the standard k-ω model, the eddy dissipation term, 𝜀, from the k-ε model is replaced with a turbulence frequency term, ω. The development of this model is described in Wilcox [55] and is noted for providing superior treatment of the viscous near-wall region as well as the effects of adverse pressure gradients. However, it is less effective in its treatment of non-turbulent free stream boundaries. The SST model, developed by Menter [56], attempts to combine the best aspects of the k-ε and k-ω models into a single formulation. This is accomplished by introducing a blending function to the ω equation, which activates the standard k-ω near the walls and activates the k-ε model away from the walls. The formulation of this 2-equation model is described in Equations 15 and 16 and the new 𝑢𝜏 is defined in Equation 17. Additionally, ANSYS CFX employs automatic wall functions which apply standard wall functions at 𝑦 + > 11 and transitions to integrating to the 37 wall at 𝑦 + < 11. Consequently, the SST model can handle any 𝑦 + value. However, to take full advantage of this model, particularly the superior near wall 𝑘-𝜔 treatment, it is recommended to aim for a 𝑦 + < 11 and at least 10 grid points in the boundary layer [53]. However, this can become computationally expensive and harder to achieve at sufficiently high Reynolds numbers common in many turbomachinery applications. The formulation of this 2-equation model and the turbulence viscosity relationship in this model is described below: 𝜕𝜌𝑘 𝜕(𝜌𝑈𝑗 𝑘) 𝜕 𝜇𝑡 𝜕𝑘 (15) + = [(𝜇 + ) ] + 𝑃 − 𝐶𝜇 𝜌𝑘𝜔 + 𝑃𝑘𝑏 𝜕𝑡 𝜕𝑥𝑗 𝜕𝑥𝑗 𝜎𝑘 𝜕𝑥𝑗 𝜕𝜌𝜔 𝜕(𝜌𝑈𝑗 𝜔) 𝜕 𝜇𝑡 𝜕𝜔 𝜔 (16) + = [(𝜇 + ) ] + 𝛼3 𝑃 + 𝜕𝑡 𝜕𝑥𝑗 𝜕𝑥𝑗 𝜎𝜔3 𝜕𝑥𝑗 𝑘 1 𝜕𝑘 𝜕𝜔 2(1 − 𝐹1 ) − β3 𝜌𝜔2 + 𝑃𝜔𝑏 𝜎𝜔2 𝜔 𝜕𝑥𝑗 𝜕𝑥𝑗 𝑘 (17) 𝑢𝜏 = 𝜌 𝜔 1 𝜕𝑘 𝜕𝜔 𝐹1 is the blending function and 𝜎 is the cross-diffusion derivative term. Further detail of 𝜔2 𝜔 𝜕𝑥𝑗 𝜕𝑥𝑗 this numerical analysis, including the development of the empirical constants, is described in [53] 3.5. Grid Generation Developing an appropriate grid is one of the most important aspects of carrying out the computational simulations. The grid must be fine enough to accurately capture the flow characteristics of the seals. However, as the mesh size increases, the computational time inevitably increases. Therefore, it is computationally advantageous to procure the minimum amount of a nodes that does not sacrifice computational accuracy. Additionally, when modelling in CFD, a proper 𝑦 + value, suitable to the chosen turbulence model is an important consideration for solution 38 accuracy. Therefore, a grid sensitivity is carried out on the gas labyrinth seals to achieve a solution that is both accurate and time efficient. For many aerodynamic applications, it is common practice to model a single sector (pie slice) of the model and apply periodic boundary conditions so that the flow through the compressor will match the simulated region. This allows for simulations to be run in a matter of minutes to a few hours on a single computer. This approach is adopted in this study when seeking mesh independents for aerodynamic parameters such as leakage and velocity profiles. The grid quality and 𝑦 + values are controlled by dividing the seal into multiple regions and refining the corresponding blocks into a sufficient number of grid elements. This also allows for better control over the transition from the boundary layer to the main flow path. Figure 15 shows the regions of mesh refinement in the axial and radial direction for three seal configurations, 8-TOR [28], 5- stepped [17], and 5-TOS, [37]. Note that Figure 15a depicts the TOS configuration. This figure also applies to the regions for the TOR configuration. Although, the SST turbulence model criteria recommends a 𝑦 + < 11, this range was only achievable for the low Reynolds number pressure condition of the 5-tooth stepped labyrinth seal. High pressure conditions correspond to high Reynolds numbers which requires a sufficiently thin first boundary layer thickness to capture the near wall boundary layer. Such high Reynolds numbers requires a drastic increase in the mesh size in order to resolve the boundary layer to a 𝑦 + < 11, while maintaining appropriate mesh statistics. However, for the higher Reynolds number cases presented in this work, there was little 𝑦 + sensitivity between 20 < 𝑦 + < 50, therefore, the labyrinth seal 𝑦 + values were kept in this range. Predicting rotordynamic performance in gas labyrinth seals requires the investigation of a 360° three-dimensional flow profiles in the compressor. This is because the pressure profile around the circumference of the seal is non-uniform. Therefore, the circumferential node count plays a 39 crucial role in determining the final grid size. As shown in Figure 16-17, the parameters used to establish the circumferential node count limit are the direct stiffness (𝐾𝑥𝑥 ) and the cross coupled stiffness (𝑘𝑥𝑦 ) . The final node count in the circumferential directions also included in Table 4. The final meshes created for the CFD simulations for three labyrinth seal configurations are shown in Figure 18. Figure 15: Regions of node refinement (a) TOS and TOR (b) Stepped labyrinth seal Table 4: Number of nodes in Figure 15 for each seal 8-TOR 5-Stepped 5-TOS 1 17 17 17 2 12 20 12 3 35 40 26 4 46 51 51 5 -- 25 -- Circumferential 160 144 144 Node Count 4110080 3693168 2922912 40 160 140 120 100 8 TOR (Kxx) KN/m 80 8 TOR (kxy) 5-Stepped (Kxx) 60 5-Stepped (kxy 40 20 0 100 125 150 175 200 225 250 Circumferential nodes Figure 16: Circumferential node study for 8 TOR and 5-Stepped labyrinth seal 100 125 150 175 200 225 250 -300 -2300 5 TOS (kxy) -340 5 TOS (Kxx) -2550 -380 kxy (KN/m) Kxx (KN/m) -2800 -420 -3050 -460 -3300 -500 -3550 -540 -3800 Circumferential nodes Figure 17: Circumferential node study for 5 TOS labyrinth seal 41 (a) (b) (c) Figure 18: Cavity mesh a) 8-TOR b) 5-Stepped c) 5-TOS 3.6. Steady State CFD The steady-state CFD method provides a well-established approach more sophisticated than bulk flow models, while maintaining an acceptable simulation turnaround. To predict the dynamic properties of the labyrinth seals, the steady state CFD method uses a whirling rotor model, which was developed by Athavale et al. [36]. It has been subsequently detailed in many other works mentioned in the literature [35, 37, 41, 42, 52]. Since every rotating shaft in a labyrinth seal rotates with an eccentricity, 𝜀, due to bending and oil-film elasticity [3], the steady state model 42 introduces an offset rotor and asymmetric fluid domain. As the rotational speed increases, the rotor undergoes a whirling motion around radius, ε, while rotating around the rotor center. In the presented CFD models, the eccentricity ε, is prescribed to be between 10-50% of the shroud clearance for each labyrinth seal. [48] notes that the seal will behave in a linear manner up to roughly 50% seal clearance. Furthermore, to reduce computational effort, the whirling rotor problem, which is inherently unsteady, is transformed into a steady-state one, by implementing a frame transfer. This process is displayed in Figure 19. Observing the motion of a rotor-seal system from a stationary frame portrays two angular velocity components, namely, the rotor spinning speed ω, and the whirling speed Ω. With the location of the rotor continually changing, the mesh, therefore, is constantly changing, which would require a transient analysis. This is avoided by transferring the stationary frame of reference into a rotating frame. In this frame, the stator wall moves in the opposite direction, or -Ω. In Ansys CFX, the stationary walls are prescribed as a counter-rotating wall. This provides a steady-state problem that permits whirl to be simulated. Moreover, the static offset in the radial direction results in a restoring force enacted upon the rotor. These forces can be determined by integrating the static pressure along the seal in the circumferential direction. 𝐹𝑥 = ∬ 𝑃𝑐𝑜𝑠(𝜃)𝑑𝐴 (18) 𝐹𝑦 = ∬ 𝑃𝑠𝑖𝑛(𝜃)𝑑𝐴 (19) The force-motion relationship (for small motion) is established using the linear system of stiffness, damping, inertia terms is mentioned in Equation 6. The x-y frame represents two orthogonal displacement directions for the of the rotating seal, relative to the stator. 43 Figure 19: Stationary to rotating frame transfer The response of the system generates an impedance curve over varying whirl frequencies. The result of this analysis is a set of impedance coefficients that quantify the forces on the labyrinth seal. A 1st order (linear) or 2nd order (quadratic) curve fit as a function of whirl speed can then be used to obtain the rotordynamic coefficients using Equations 20 and 21. In this frame, the radial impedance force is in the same direction as 𝐹𝑥 , and the tangential impedance force is in the same direction as 𝐹𝑦 . 𝐹𝑟 𝐹𝑥 (20) = = −𝐾𝑥𝑥 − 𝑐𝑥𝑦 Ω + 𝑀𝑥𝑥 Ω2 𝜀 𝜀 𝐹𝑡 𝐹𝑦 (21) = = 𝑘𝑥𝑦 − 𝐶𝑥𝑥 Ω − 𝑚𝑥𝑦 Ω2 𝜀 𝜀 3.7. Transient CFD In an alternative approach, a time-marching transient whirling rotor model can be implemented to obtain rotordynamic coefficients. Unlike the steady-state approach, this method does not necessitate an axisymmetric geometry, or a frequency-independent assumption, but can be extended to non-axisymmetric geometry and frequency-dependent rotordynamic characteristics. Since there is no coordinate frame transformation in this approach, the seal is 44 initially maintained concentric, and is physically perturbed into a periodic whirling rotor orbit using moving mesh techniques. CFX does not allow the coordinate frame to move, therefore, this requires the velocities on the rotor surface to be defined in the global coordinate frame. To illustrate this, Figure 20 is used to derive the rotor velocities. From Figure 20, the x and y velocity components are derived in Equations 22-27. The negative sign for 𝑉𝑥 corresponds to a right-handed coordinate system definition. 𝑉𝑥 = −𝑟𝜔sin (𝜃) (22) 𝑉𝑦 = 𝑟𝜔cos (𝜃) (23) And, 𝑥 − 𝑥𝑟 (24) cos(𝜃) = 𝑟 𝑦 − 𝑦𝑟 (25) 𝑠𝑖𝑛(𝜃) = 𝑟 Therefore, the final expression for the surface velocities can be written as: 𝑉𝑥 = −𝜔(𝑦 − 𝑦𝑟 ) (26) 𝑉𝑦 = 𝜔(𝑥 − 𝑥𝑟 ) (27) The restoring forces are then calculated from the fluid acting upon the rotor surface and the impedances can be determined from the linear force-motion relationship mentioned in Equation 6. Moreover, the extraction of rotordynamic coefficients using the transient approach is dependent up the whirling rotor model implemented. These whirling orbits models are mentioned in the following section. 45 Figure 20: Derivation of the Rotor Surface Velocity 3.7.1. Rotor Whirling Models This research considers whirling models with two degrees of freedom (2DOF), in which the rotor moves in two orthogonal directions, x, and y. First, a circular whirling rotor model imposes a periodic circular orbit in a transversal plane normal to the rotational axis. This process is depicted in Figure 21. A transient simulation can be run for a single excitation frequency or multiple excitation frequencies. For a single frequency whirling motion, the displacements are defined in Equations 28-29 [47]. Note that Equations 28-29 describe a forward whirling rotor. Equations 30-31 applies a backward whirling rotor model. Since a single frequency whirling model contains only one whirling frequency in each simulation, it requires multiple transient simulations for each whirl frequency. To allow for direct extraction of rotordynamic coefficients requires a transient simulation for forward and backward whirl at the same frequency. 46 Forward whirl 𝑥 = 𝜀 ∗ 𝑐𝑜𝑠 (Ω𝑡) (28) 𝑦 = 𝜀 ∗ 𝑠𝑖𝑛 (Ω𝑡) (29) Backward whirl 𝑥 = 𝜀 ∗ 𝑐𝑜𝑠 (Ω𝑡) (30) 𝑦 = −𝜀 ∗ 𝑠𝑖𝑛 (Ω𝑡) (31) A multi-frequency whirling model is an attempt to reduce the computational time by introducing multiple whirl frequencies in the periodic whirling rotor, by applying a linear combination of Equations 32-35. One multi-frequency transient simulation may take longer than the single frequency simulation, however, only two transient simulations are needed to extract rotordynamic coefficients over a range of excitation frequencies. In this approach, the circular whirling orbit is expressed as: Forward whirl 𝑁 (32) 𝑥 = 𝜀 ∑ 𝑐𝑜𝑠 (Ω𝑘 𝑡) 𝑘=1 𝑁 (33) 𝑦 = 𝜀 ∑ 𝑠𝑖𝑛 (Ω𝑘 𝑡) 𝑘=1 Backward whirl 𝑁 (34) 𝑥 = 𝜀 ∑ 𝑐𝑜𝑠 (Ω𝑘 𝑡) 𝑘=1 𝑁 (35) 𝑦 = −𝜀 ∑ 𝑠𝑖𝑛 (Ω𝑘 𝑡) 𝑘=1 47 Figure 21: Circular whirl orbit a) forward whirl, b) backward whirl An additional method applied in this work uses the same transient technique; however, a periodic elliptical whirling motion is applied instead of a circular orbit. In this case, the rotor whirls about a major and minor axis, a and b respectively, on the coordinate axis as depicted in Figure 22. The major axis represents the direction in which the rotor is excited. Therefore, expressions for single and multi-frequency elliptical orbit motion can be written in the following form: Single-frequency: x-direction excitation 𝑥 = 𝑎𝑐𝑜𝑠 (Ω𝑡) (36) 𝑦 = 𝑏𝑠𝑖𝑛 (Ω𝑡) (37) Single-frequency: y-direction excitation 𝑥 = 𝑏𝑐𝑜𝑠 (Ω𝑡) (38) 𝑦 = 𝑎𝑠𝑖𝑛 (Ω𝑡) (39) 48 Multi-frequency: x-direction excitation 𝑁 (40) 𝑥 = 𝑎 ∑ 𝑐𝑜𝑠 (Ω𝑘 𝑡) 𝑘=1 𝑁 (41) 𝑦 = 𝑏 ∑ 𝑠𝑖𝑛 (Ω𝑘 𝑡) 𝑖=1 Multi-frequency: x-direction excitation 𝑁 (42) 𝑥 = 𝑏 ∑ 𝑐𝑜𝑠 (Ω𝑘 𝑡) 𝑘=1 𝑁 (43) 𝑦 = 𝑎 ∑ 𝑠𝑖𝑛 (Ω𝑘 𝑡) 𝑘=1 Figure 22: Elliptical whirl orbit a) x-direction excitation, b) y-direction excitation The whirl frequency range Ω𝑘 begins at the prescribed minimum frequency, described as the base frequency, and N is the frequency number which depends on the frequency range implemented in the CFD simulation. Ω𝑘 can be expressed as Ω𝑘 = 2𝜋𝑓𝑘 , where 𝑓𝑘 is the rotor vibration frequency. 49 The base frequency can be an estimation of the cross-over frequency or the damped natural frequency of the rotor [57]. The maximum frequency, and thus, the frequency range is related to the time-step implemented in the CFD. To apply the FFT, the number of timesteps in the base frequency must be a power of 2. Therefore, the given minimum and maximum frequency is related to the time step in the following Equations, which are adapted from [57]. 1 (44) ∆𝑡 ≤ 𝜏Ω𝑚𝑎𝑥 1 (45) = 2𝑛 ∆𝑡Ω𝑚𝑖𝑛 Where 𝜏 is the number of timesteps per revolution which is the chosen number of timesteps in a single period of the maximum frequency component. The choice of the timestep on the transient solution will be further discussed as this work progresses. 3.7.2. Transient Rotordynamic Coefficient Identification . Two different approaches are used to post-process the data obtained from the transient CFD data. First, rotordynamic coefficients can be extracted by implementing a KCM model fit which assumes frequency independence, where the mass, stiffness, and damping coefficients remain constant for a given running speed. In this model, a quadratic curve is fit to the real and imaginary parts of the force divided by the displacement [50]. For a periodic circular whirling orbit, the orthogonal force components, 𝐹𝑥 and 𝐹𝑦 , are transferred to the radial and tangential force in Equations. 𝐹𝑟 = 𝐹𝑥 𝑐𝑜𝑠(Ω𝑡) + 𝐹𝑦 sin (Ω𝑡) (46) 𝐹𝑡 = −𝐹𝑥 𝑠𝑖𝑛(Ω𝑡) + 𝐹𝑦 cos (Ω𝑡) (47) And substituting Equations 46-47 into the linear force-motion Equation 3, 𝐹𝑥 and 𝐹𝑦 can be substituted with the stiffness, mass, and damping terms. Inspection of Equations 46-47 shows that 50 when rotor whirls to a position Ω𝑡 = 0° , 𝐹𝑟 = 𝐹𝑥 and 𝐹𝑡 = 𝐹𝑦 . Moreover, when the rotor whirls to a position of Ω𝑡 = 90° , 𝐹𝑟 = 𝐹𝑦 and 𝐹𝑡 = −𝐹𝑥 . Therefore, the rotordynamic coefficients can be extracted from a linear or parabolic curve fit from impedance plots for forward whirl using Equations 48-51. Moreover, at least three simulations at three different orbital frequencies are needed to obtain a KCM curve fit for the rotordynamic coefficients. 𝐹𝑥,0 (48) = −𝐾𝑥𝑥 − 𝑐𝑥𝑦 Ω + 𝑀𝑥𝑥 Ω2 𝜀 𝐹𝑥,90 (49) = 𝑘𝑥𝑦 − 𝐶𝑥𝑥 Ω − 𝑚𝑥𝑦 Ω2 𝜀 𝐹𝑦,0 (50) = −𝑘𝑦𝑥 − 𝐶𝑦𝑦 Ω + 𝑚𝑦𝑥 Ω2 𝜀 𝐹𝑦,90 (51) = −𝐾𝑦𝑦 + 𝐶𝑦𝑥 Ω + 𝑀𝑦𝑦 Ω2 𝜀 A similar operation is carried out for an elliptical whirling rotor model. However, taking the derivative and second derivative of Equations and substituting into the linear force-motion relationship gives a set of equations. Given the vibration amplitudes in the major and minor axis, a set of equations at Ω𝑡 = 0° and Ω𝑡 = 90° , 𝐹𝑥,0 = −𝑎𝐾𝑥𝑥 − 𝑏𝑐𝑥𝑦 Ω + 𝑎𝑀𝑥𝑥 Ω2 (52) 𝐹𝑥,90 = 𝑏𝑘𝑥𝑦 − 𝑎𝐶𝑥𝑥 Ω − 𝑏𝑚𝑥𝑦 Ω2 (53) 𝐹𝑦,0 = −𝑎𝑘𝑦𝑥 − 𝑏𝐶𝑦𝑦 Ω + 𝑎𝑚𝑦𝑥 Ω2 (54) 𝐹𝑦,90 = −𝑏𝐾𝑦𝑦 + 𝑎𝐶𝑦𝑥 Ω + 𝑏𝑀𝑦𝑦 Ω2 (55) 51 An alternative approach is to transfer the time-domain into a frequency domain using a Fast Fourier Transform (FFT). To accomplish this, Equation 6 is transformed into the following frequency domain in the following matrix notation. 𝐹𝑥 𝐻𝑥𝑥 𝐻𝑥𝑦 𝐷𝑥 (56) [𝐹 ] = [ ][ ] 𝑦 𝐻𝑦𝑥 𝐻𝑦𝑦 𝐷𝑦 𝐹𝑥 and 𝐹𝑦 are the fluid response forces acting upon the rotor and 𝐷𝑥 and 𝐷𝑦 are the frequency domain components of the relative rotor displacement in the x and y direction. Both of these components are obtained directly from monitoring the transient mesh motion and force responses in the CFD solver. Equation gives 2 equations with 4 unknowns. Therefore, two separate transient runs are necessary for the 2DOF whirling rotor models in order to obtain four equations for the four unknowns. For the circular whirling rotor model this includes forward and backward whirling excitation. For the elliptical whirling model, this includes excitation in the x-direction and excitation in the y-direction. 𝐹𝑥𝑥 𝐹𝑥𝑦 𝐻𝑥𝑥 𝐻𝑥𝑦 𝐷𝑥𝑥 𝐷𝑥𝑦 (57) [ ]=[ ][ ] 𝐹𝑦𝑥 𝐹𝑦𝑦 𝐻𝑦𝑥 𝐻𝑦𝑦 𝐷𝑦𝑥 𝐷𝑦𝑦 With this information, the force impedances in the frequency domain can be obtained, which provide the frequency-dependent stiffness and damping coefficients. 𝐻𝑖𝑗 = 𝐾𝑖𝑗 + 𝑗(𝐶𝑖𝑗 Ω) − 𝑀𝑖𝑗 Ω2 (58) Where 𝑗 = √−1. Therefore, this equation can be divided into a real and imaginary part as follows. 𝑅𝑒(𝐻𝑖𝑗 ) = 𝐾𝑖𝑗 − 𝑀𝑖𝑗 Ω2 (59) 𝐼𝑚(𝐻𝑖𝑗 ) = 𝑗(𝐶𝑖𝑗 Ω) (60) For gas labyrinth seals 𝑀𝑖𝑗 is small and is usually neglected. It is important to note that for the force response and mesh displacement monitored in CFX, 𝐹𝑖𝑗 and 𝐷𝑖𝑗 , the first subscript denotes the whirling direction, and the second subscript denotes the direction of the force response 52 direction. For the rotordynamic coefficients, 𝐾𝑖𝑗 and 𝐶𝑖𝑗 , the first subscript denotes the direction of the force response, and the second subscript denotes the whirling direction. Finally, the frequency-dependent stiffness and damping coefficients can be directly extracted from the real and imaginary components. 𝐾𝑖𝑗 = 𝑅𝑒(𝐻𝑖𝑗 ) (61) 𝐼𝑚(𝐻𝑖𝑛 ) (62) 𝐶𝑖𝑗 = Ω 53 Chapter 4: Steady State CFD Results 4.1. 8-Tooth on Rotor Labyrinth Seal The numerical results for the steady state analysis of the 8-tooth on rotor labyrinth seal are presented here. Figures 23 and 24 show the predicted radial and tangential impedances for the 8- tooth on rotor labyrinth seal. The impedance forces are obtained for two different inlet swirl values. Additionally, the k-ε and SST turbulence models are both used and compared to the experimental data. A total of 5 whirl speeds are solved for the rotordynamic coefficients. The tangential impedance curve produces essentially linear results. Therefore, the linear experimental curve (no inertial coefficient) is included in the graph for comparison. However, the CFD produces a radial impedance curvature that fits well with a second-order curve fit. This is contrary to the experiment, which details a linear relationship. This difference can be explained by the small range of excitation frequencies used in Pelletti [28], which was too small to capture a 2nd order curve. Therefore, a linear experimental curve fit is not included in the radial impedance plot. However, the added inertial term in the predicted radial impedance is small. Moreover, a quantitative comparison of the rotordynamic coefficients derived from each curve fit at both inlet swirl ratios is presented in Tables 5-6. For the no swirl condition, the predicted rotordynamic coefficients provide close agreement with the rotordynamic coefficients measured in experiment. The steady-state models are comparable in direct stiffness, which is slightly overpredicted, and direct damping, which is somewhat underpredicted. At this condition, there is little advantage of using one turbulent model over the other, although the k-ε turbulence model gives a marginally better prediction of the cross- coupled stiffness terms. However, at the higher inlet swirl condition, the difference between the turbulence models becomes more pronounced. As Table 6 shows, the SST model achieves a better agreement with the experimental direct and cross-coupled stiffness value, whereas the k-ε model 54 matches slightly better with the direct damping term. For the SST model, 𝐶𝑥𝑥 increases with swirl, whereas the experiment finds a decrease in 𝐶𝑥𝑥 as swirl increases. From these observations, it can be concluded that the choice of turbulence modelling becomes more important when considering higher inlet swirl ratios. Figure 23: Radial impedance for the TOR labyrinth at 0 inlet swirl 55 Figure 24: Tangential impedance for the TOR labyrinth at 0 inlet swirl Figure 25: Radial impedance for the TOR labyrinth at 0.165 inlet swirl 56 Figure 26: Tangential impedance for the TOR labyrinth at 0.165 inlet swirl Table 5: Rotordynamic coefficients For TOR labyrinth seal at 0 inlet swirl ratio k-ε model SST model Experiment 𝐾𝑥𝑥 (N/m) -118250 -116571 -146910 𝑐𝑥𝑦 (N-s/m) 100 82 70 𝐶𝑥𝑥 (N-s/m) 395 397 453 𝑘𝑥𝑦 (N/m) 107798 129551 112542 𝑀𝑥𝑥 (kg) 0.073 0.059 𝑚𝑥𝑦 (kg) 0.0087 0.045 Table 6: Rotordynamic coefficients For TOR labyrinth seal at 0.165 inlet swirl ratio k-ε model SST model Experiment 𝐾𝑥𝑥 (N/m) -128,977 -116103 -88023 𝑐𝑥𝑦 (N-s/m) 111 80 -104 𝐶𝑥𝑥 (N-s/m) 341 410 202 𝑘𝑥𝑦 (N/m) 179660 231121 323663 𝑀𝑥𝑥 (kg) 0.083 0.059 𝑚𝑥𝑦 (kg) 0.036 0.03 57 4.2. 5-Stepped Labyrinth Seal The analysis presented here applies to the stepped labyrinth seal configuration characterized by 5 teeth on the stationary domain. The rotordynamic coefficients reported in Kwanka [17] were measured in relation to the circumferential velocity entering the seal. Figures 27a-d show the results for rotordynamic stiffness and damping coefficients, which are plotted versus circumferential velocity Cu. Gao [52] also performed a CFD analysis on the same stepped labyrinth seal. The results published in this work are also included in the graphs for comparison. Moreover, the effects of turbulence modeling were also considered, comparing the rotordynamic coefficients using k-ε and k-ω based SST turbulence models. Figure 27a shows that the CFD underpredicts the direct stiffness values compared to experimental measurements. In most cases, the CFD direct stiffness values decrease as the swirl velocity increases. This is found to be true for the presented k-ε model and the CFD model reported by Gao [41]. However, this trend is not as clear for the SST model as the direct stiffness begins to increase at significantly high swirl velocities. Contrary to the CFD, the direct stiffness values from experiment were positive and increased as the inlet circumferential velocity increase. Overall, the present k-ε model utilized in this analysis provided the closest realization of direct stiffness values with experimental data. The calculated cross-coupled stiffness values presented in Figure 27b. correctly predict an increase in the 𝑘𝑥𝑦 value as the inlet circumferential velocity increases. In Gao’s model, the 𝑘𝑥𝑦 values matched close to experiment for lower circumferential velocity conditions. This difference increased dramatically as the inlet swirl increased. In comparison to Gao [52], both the k-ε and SST models predicted slightly higher cross-coupled stiffness at low circumferential velocity conditions but came into a better agreement at higher inlet swirl conditions. The SST model predicted slightly closer cross-coupled stiffness than the k-ε model. It 58 is seen in Figure 27c-d that damping terms. Moreover, the present CFD models in this study shows improvement in 𝐶𝑥𝑥 and 𝑐𝑥𝑦 when compared to [52]. the inlet circumferential velocity has little influence on the direct damping and cross-coupled. Figure 27: Comparison of rotordynamic coefficients for 5-tooth stepped seal (a) direct stiffness (b) cross-coupled stiffness (c) direct damping (d) cross-coupled damping 4.3 5-Tooth on Stator Labyrinth Seal To establish the consistency of steady-state CFD numerical models, a 3D 5 TOS seal-only model is compared with results from existing CFD and bulk flow models [37]. For this purpose, the k-ε turbulence model with scalable wall functions was chosen for this analysis to match operating conditions with previous work. A no pre-swirl condition was considered. The radial and tangential impedance impedances for this case are displayed in Figures 28-29. Figure 29 shows 59 that for the tangential impedance curve, the CFD provides consistent results with the CFD models reported in the literature. Moreover, the bulk flow model prediction also provides close validation. The radial impedance curve in Figure 28 shows the closest agreement with Hirano et al. [37]. The discrepancy between TASCFlow model can potentially be attributed to the difference in grid resolution. A courser mesh was used in [37], which had 761,852 nodes compared to 2,922,912 nodes used in the present model. The CFD results show a larger discrepancy with the bulk flow program and the CFD results reported in [41, 42]. Gao [41] and Subramanian et al. [42] both concluded that the difference in the radial impedance values could be neglected because the direct stiffness and cross-coupled damping terms have a negligible influence on rotor stability. This is generally true for throughflow compressors [18]. However, large stiffness values can indirectly affect rotor stability by influencing the vibrational frequencies that affect the rotor's critical speed. This has been reported for longer labyrinth seals [18, 28]. In this situation, accurately predicting 𝐾𝑥𝑥 becomes important. Therefore, further investigation with experimental measurements is needed to shed light on the discrepancy found between the steady-state CFD and bulk flow models. An additional option for further investigation would be to use unsteady CFD. A quantitative comparison of the 1st order rotordynamic coefficients for each method is presented in Table 7. 60 Figure 28: Radial impedances for 5-TOS labyrinth at 0 inlet swirl Figure 29:Tangential impedances for 5-TOS labyrinth at 0 inlet swirl Table 7: Comparison of predicted seal coefficients with the results reported by [37, 41, 42] 𝐾𝑥𝑥 𝑐𝑥𝑦 𝐶𝑥𝑥 𝑘𝑥𝑦 CFD (Present) 2815069 978 1927 -524945 Hirano CFD [37] 1880000 1240 1600 -720000 Gao CFD [41] 6010000 1330 1890 -856000 Subramanian CFD [42] 650000 700 1300 -580000 Bulk flow [37] 375000 567 1050 -400000 61 4.4. Conclusions Up to this point, the steady state CFD method showed good agreement with available experimental data. For the 8-TOR labyrinth seal, the presented k-ε and SST turbulence models gave comparable results for the no swirl case and showed close agreement to experiment. For the higher swirl case, the effect of turbulence model was more noticeable. However, additional case studies are needed to determine the superiority of chosen turbulence models. For the 5-teeth stepped seal, the CFD failed to predict the correct sign of the direct stiffness for all but the 0-swirl case in the presented k-ε model. At sufficiently high inlet swirl, the SST model provides better prediction for the cross-coupled stiffness term. For the 5-TOS labyrinth seal, the CFD models presented by each author provided comparable results for the tangential impedance forces. However, a larger discrepancy is observed for the predicted radial impedance force, giving a pessimistic comparison amongst the authors. Therefore, further investigations will be required to determine the inconsistency between the authors when predicting the radial impedances. Transient CFD may shed light on this area. 62 Chapter 5: Steady State vs Transient CFD Comparison In the past, rotordynamic analysis on labyrinth seals has been dominated by simple in- house bulk flow models, and steady state CFD analysis. The previous chapter demonstrated the accuracy of using these techniques. However, CFD has now matured to the point where the time required to perform unsteady analyses on these components is becoming realistic in the common work setting. Adding the time-marching variables inherent in transient solutions gives a potential accuracy advantage. This chapter will comprehensively compare steady state and the transient CFD techniques used to predict the aerodynamic performance and rotordynamic coefficients in gas compressor labyrinth seals. Part of this goal weighs the benefits between computational time and accuracy. To realize this goal, a full 3D eccentric steady state model described in chapter 3.6 and two 3D whirling rotor transient models, mentioned in chapter 3.7 are implemented in this work. To accomplish this effort, experimental data provided from the literature is used for comparison. Furthermore, the influence of pre-swirl is also considered. 5.1. CFD Model and Set up A recent rotordynamic experiment carried out at Texas A&M provides rotordynamic coefficients for a 14 TOR labyrinth seal at high pressure conditions. This section expands upon the modelling set up for this seal. The turbulent model chosen for this study is the 𝜅-𝜀 model. The transient and steady CFD models are first verified by comparing the predicted leakage value with the experimental leakage in Figure 31. Experimental leakage values were calculated from measured volumetric flow, temperatures, and pressures referenced in equation 8 and 9 of Arthur [22]. Furthermore, this comparison is analyzed at different pre-swirl ratios at 10,200 rpm. Figure 31 demonstrates that each CFD method gives identical leakage predictions. This preliminary 63 0.2 Steady State CFD 0.19 Transient CFD Leakage [kg/s] Experiment 0.18 0.17 0.16 0.15 -0.01 0.09 0.19 0.29 0.39 0.49 0.59 Preswirl Ratio Figure 30: Comparison of Experimental Leakage measurements with steady and transient CFD investigation demonstrates that a faster steady state CFD analysis is sufficient when leakage is the only calculation of interest. Additionally, the experiment shows an increase in leakage towards medium inlet swirl, and then decreasing at high inlet swirl. This trend is not typical. The CFD does not capture this trend, but rather shows a decrease in leakage as swirl increases, which is consistent with trends found in previous literature. Nevertheless, this difference is small and both CFD methods are in close agreement with the experimental leakage values. When modelling a labyrinth seal-only model, it is often necessary to implement artificial inlets and outlets to achieve convergence within the solver. This is because the flow recirculates just upstream the first tooth of the labyrinth seal. Extending the inlet and outlet allows for the flow field to resolve and converge. However, this can add a considerable length upstream the first tooth. Therefore, the swirl entering at the inlet must travel an extended distance before entering the first tooth. To compensate for this, a higher swirl ratio is implemented at the extended inlet to achieve a matching swirl ratio at the location entering the seal recorded in the experiment. Figure 32 shows 64 the average swirl ratio traveling through the labyrinth seal. The effects of inlet swirl last up to the midpoint of the labyrinth seal before becoming uniform for all inlet swirl conditions. 0.8 Low 0.7 Preswirl 0.6 Medium Swirl 0.5 Swirl Ratio High 0.4 Preswirl 0.3 0.2 0.1 0 -4 6 16 26 36 46 56 66 76 Axial Location [mm] Figure 31: Swirl ratio through each seal cavity 5.2. Post Processing This section describes the CFD models used to predict the rotordynamic coefficients for the 14 TOR labyrinth seal. The whirling rotor simulation is modeled in two ways. The first method is mentioned in section 3.6. which uses the steady state CFD approach, transforming the stationary frame of reference into a rotating frame. To simulate whirl, the geometry of the rotor is offset by a chosen eccentricity. For the 14-tooth on rotor labyrinth seal, an eccentricity of 0.02 mm which corresponds to 20% of the labyrinth seal clearance is chosen for the steady state model. This value is chosen to capture the linear, small motion characteristics of the labyrinth seal. After the simulation setup, the model is run at a series of different whirl frequencies defined as 1x, 0.75x, 0.50x, and 0.25x of the running speed. The whirl frequencies were chosen to allow for an appropriate number of points to create the model curve fits. Additionally, a static eccentric case is also simulated, which offsets the rotor at the prescribed eccentricity without whirl. 65 The second CFD method uses a time-marching whirling rotor model. Since there is no coordinate frame transformation in this approach, the seal is initially concentric. This simulation imposes a whirling rotor wall and then calculates the restoring forces from the fluid on the rotor surface. In CFD, this requires a moving mesh. The results obtained from the CFD data are post-process using the two different approaches mentioned in chapter 3.7. First, rotordynamic coefficients are extracted by implementing a KCM model fit which assumes frequency independence. In this model, a quadratic curve is fit to the real and imaginary parts of the force divided by the displacement [50]. The same whirling frequencies as the steady state model are chosen for this model, and only a single whirling frequency is run for each simulation. Transient parameters 2𝜋 𝑟𝑎𝑑 selected for this study include a time step of 𝑟𝑎𝑑 = 6.536 × 105 with 4 iterations per 90∗(1062 ) 𝑠 timestep. This corresponds to 90 timesteps per revolution. To determine convergence for transient whirl cases, the force profile is monitored until it becomes periodic. A periodic solution is achieved after about 3-4 revolutions. This means that the force difference between the nearby periodicities is less than 0.2%, which is suggested by [45]. For the static 0-whirl case, convergence is achieved when the force profile flattens out. Furthermore, the transient whirling simulations is started from the steady-state concentric simulations. By running a concentric steady state simulation first, it is verified that for zero eccentricity, the net forces acting upon the rotor are zero. The second transient approach used to compute rotordynamic coefficients generates data for two linearly independent excitations and solves for the matrix Equation 57. In this approach, six excitation frequencies chosen from the experimental data presented by Arthur [22] for each pre-swirl case. The timestep for the multi-frequency whirling orbit model is slightly adjusted so that the standard Fast Fourier Transform algorithm can be applied. This is done by ensuring that the number of time steps in the base frequency is a multiple of 2. 66 1 (63) Δ𝑡 = 𝑓𝑜 ∗ 2𝑛 For the 14 TOR labyrinth seal, 512-time steps (2𝑛 = 512) of the base frequency (30 Hz), were required to reach a converged solution. This resulted in a time step of 6.51 × 105 s. Additionally, this requires a minimum of 1024 time steps total or 2 cycles of the base frequency, to allow for the initial startup transients to subside. Nevertheless, it is recommended that 3 cycles of the base frequency are needed to ensure periodic repeatability between cycles 2 and 3. Figures 33-69 display the monitored displacement and force responses used to determine the rotordynamic coefficients in the following section. For the single whirling frequency cases, only the periodic time domain is represented for the circle orbit. For the multi-frequency approach, the time domain is plotted along with its frequency domain representation, which depict the peak magnitudes. Figures 33-36 give examples of the monitored displacement data for circular and elliptical orbit rotor whirling motions. These graphs are the same for each inlet swirl case. The frequency domain displacement representatives are seen display in Figures 37-39. It is noticeable that the monitored displacements slightly underpredict the target displacement amplitudes. This variation is due to the numerical error associated with the mesh deformation in the transient CFD solution. However, this deviation is negligible. Therefore, for each whirling motion, the monitored displacement values can accurately capture the target vibration amplitudes predefined in Equations 32-35 & 40-43. Additionally, the whirling motion and force response figures show that the fluid response forces have the same frequency components but different initial phases with the rotor whirling motions. In the frequency domain representation, it is noticeable that the force amplitude is frequency dependent, even though the same whirling amplitude is prescribed at each frequency. At low inlet swirl ratios, the force response magnitudes increase with increasing frequency for 67 each whirling orbit (circle or elliptical). However, at higher inlet swirl ratios the force response may decrease and then increase as the frequency increases, depending about the prescribed whirling motion. For example, this occurs at medium (0.25) and high swirl (0.48) for a forward whirling motion within a circular orbit. A circular orbit with backward whirling motion demonstrates and increasing force amplitude as the excitation frequency increases. When prescribing an elliptical whirling orbit, there is little change in the force response when imposing a vibration amplitude with x excitation or y excitation. To summarize, the rotor whirling parameters such as whirling orbit, vibration amplitude, and frequency component have significant influence on the rotor whirling motions and force responses. However, as will be demonstrated in chapter 5.3., the whirling parameters have little influence on the computed rotordynamic coefficients. 8.0 x y 6.0 4.0 Displacement (𝜇m) 2.0 0.0 -2.0 -4.0 -6.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 32: Monitored whirling motion for circular orbit, forward excitation in the time domain 68 8.0 x 6.0 y 4.0 Displacement (𝜇m) 2.0 0.0 -2.0 -4.0 -6.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 33: Monitored whirling motion for circular orbit, backward excitation in the time domain 7.0 x 6.0 y 5.0 Displacement (𝜇m) 4.0 3.0 2.0 1.0 0.0 -1.0 -2.0 -3.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 34:Monitored whirling motion for elliptical orbit, x excitation in the time domain 69 5.0 x y 4.0 3.0 2.0 Displacement (𝜇m) 1.0 0.0 -1.0 -2.0 -3.0 -4.0 -5.0 0.00 0.02 0.04 0.06 0.08 0.10 Time Step Figure 35: Monitored whirling motion for elliptical orbit, y excitation in the time domain 1.20 x y Target Amplitude 1.00 0.80 Displacement (𝜇m) 0.60 0.40 0.20 0.00 0 30 60 90 120 150 180 Frequency (Hz) Figure 36: FFT of whirling motion for circular orbit, forward and backward excitation 70 x 1.20 y Target Amplitude 1.00 0.80 Displacement (𝜇m) 0.60 0.40 0.20 0.00 0 30 60 90 120 150 180 Frequency (Hz) Figure 37: FFT of whirling motion, elliptical orbit, x excitation x 1.20 y Target Amplitude 1.00 0.80 Displacement (𝜇m) 0.60 0.40 0.20 0.00 0 30 60 90 120 150 180 Frequency (Hz) Figure 38: FFT of whirling motion, elliptical orbit, y excitation 71 125 Fx Fy 75 Reaction Force (N) 25 -25 -75 -125 0.000 0.010 0.020 0.030 0.040 Time (s) Figure 39: Monitored force response for single-frequency circular orbit, backward excitation. 0.08 preswirl. 40 Fx Fy 30 20 Reaction Force (N) 10 0 -10 -20 -30 -40 0.00 0.01 0.02 0.03 0.04 Time (s) Figure 40: Monitored force response for single-frequency circular orbit, forward excitation 0.08 preswirl. 72 15.0 Fx 10.0 Fy 5.0 Reaction Force (N) 0.0 -5.0 -10.0 -15.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 41: Monitored force response for multi-frequency circular orbit, backward excitation. 0.08 preswirl. 3.5 Fx Fy 3 2.5 Reaction Force (N) 2 1.5 1 0.5 0 -10 20 50 80 110 140 170 200 Frequency (Hz) Figure 42: FFT of fluid force response for circular orbit, backward excitation. 0.08 preswirl 73 Fx 15.0 Fy 10.0 5.0 Reaction Force (N) 0.0 -5.0 -10.0 -15.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 43: Monitored force response for multi-frequency circular orbit, forward excitation. 0.08 preswirl. Fx 0.9 Fy 0.8 0.7 0.6 Reaction Force (N) 0.5 0.4 0.3 0.2 0.1 0 0 30 60 90 120 150 180 Frequency (Hz) Figure 44: FFT of fluid force response for circular orbit, forward excitation. 0.08 preswirl 74 5.0 Fx 4.0 Fy 3.0 2.0 Reaction Force (N) 1.0 0.0 -1.0 -2.0 -3.0 -4.0 -5.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 45: Monitored force response for multi-frequency elliptical orbit, x excitation. 1.4 Fx 1.2 Fy 1 Reaction Force (N) 0.8 0.6 0.4 0.2 0 0 30 60 90 120 150 180 Frequency (Hz) Figure 46: FFT of fluid force response for elliptical, x excitation. 0.08 preswirl 75 5.0 Fx 4.0 Fy 3.0 2.0 Reaction Force (N) 1.0 0.0 -1.0 -2.0 -3.0 -4.0 -5.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 47: Monitored force response for multi-frequency elliptical orbit, y excitation Fx 1.4 Fy 1.2 1 Reaction Force (N) 0.8 0.6 0.4 0.2 0 0 30 60 90 120 150 180 Frequency (Hz) Figure 48: FFT of fluid force response for elliptical, y excitation. 0.08 preswirl 76 20 Fx Fy 15 10 Reaction Force (N) 5 0 -5 -10 -15 -20 0.00 0.01 0.02 0.03 0.04 Time (s) Figure 49: Monitored force response for single-frequency circular orbit, backward excitation. 0.25 preswirl. 40 Fx Fy 30 20 Reaction Force (N) 10 0 -10 -20 -30 -40 0.00 0.01 0.02 0.03 0.04 Time (s) Figure 50: Monitored force response for single-frequency circular orbit, forward excitation. 0.25 preswirl. 77 Fx 20.0 Fy 15.0 10.0 Reaction Force (N) 5.0 0.0 -5.0 -10.0 -15.0 -20.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 51: Monitored force response for multi-frequency circular orbit, backward excitation. 0.25 preswirl. Fx 4 Fy 3.5 3 Reaction Force (N) 2.5 2 1.5 1 0.5 0 0 30 60 90 120 150 180 Frequency (Hz) Figure 52: FFT of fluid force response for circular orbit, backward excitation. 0.25 preswirl 78 5.0 Fx 4.0 Fy 3.0 2.0 Reaction Force (N) 1.0 0.0 -1.0 -2.0 -3.0 -4.0 -5.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 53: Monitored force response for multi-frequency circular orbit, forward excitation. 0.25 preswirl. 0.9 Fx Fy 0.8 0.7 0.6 Reaction Force (N) 0.5 0.4 0.3 0.2 0.1 0 0 30 60 90 120 150 180 Frequency (Hz) Figure 54: FFT of fluid force response for circular orbit, forward excitation. 0.25 preswirl 79 10.0 Fx 8.0 Fy 6.0 4.0 Reaction Force (N) 2.0 0.0 -2.0 -4.0 -6.0 -8.0 -10.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 55: Monitored force response for multi-frequency elliptical orbit, x excitation. 0.25 preswirl. 1.6 Fx Fy 1.4 1.2 Reaction Force (N) 1 0.8 0.6 0.4 0.2 0 0 30 60 90 120 150 180 Frequency (Hz) Figure 56: FFT of fluid force response for elliptical orbit, x excitation. 0.25 preswirl 80 10.0 Fx 8.0 Fy 6.0 4.0 Reaction Force (N) 2.0 0.0 -2.0 -4.0 -6.0 -8.0 -10.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 57: Monitored force response for multi-frequency elliptical orbit, y excitation. 0.25 preswirl. 1.6 Fx Fy 1.4 1.2 Reaction Force (N) 1 0.8 0.6 0.4 0.2 0 0 30 60 90 120 150 180 Frequency (Hz) Figure 58: FFT of fluid force response for elliptical orbit, y excitation. 0.25 preswirl 81 150 Fx Fy 100 50 Reaction Force (N) 0 -50 -100 -150 0.00 0.01 0.02 0.03 0.04 Time (s) Figure 59: Monitored force response for single-frequency circular orbit, backward excitation. 0.48 preswirl. 30 Fx Fy 20 10 Reaction Force (N) 0 -10 -20 -30 0.00 0.01 0.02 0.03 0.04 Time (s) Figure 60: Monitored force response for single-frequency circular orbit, forward excitation. 0.48 preswirl 82 40 Fx 30 Fy 20 Reaction Force (N) 10 0 -10 -20 -30 -40 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 61: Monitored force response for multi-frequency circular orbit, backward excitation. 0.48 preswirl. 4.5 Fx Fy 4 3.5 3 Reaction Force (N) 2.5 2 1.5 1 0.5 0 0 30 60 90 120 150 180 Frequency (Hz) Figure 62: FFT of fluid force response for circular orbit, backward excitation. 0.48 preswirl 83 4.0 Fx 3.0 Fy 2.0 Reaction Force (N) 1.0 0.0 -1.0 -2.0 -3.0 -4.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 63: Monitored force response for multi-frequency circular orbit, forward excitation. 0.48 preswirl. 0.8 Fx Fy 0.7 0.6 Reaction Force (N) 0.5 0.4 0.3 0.2 0.1 0 0 30 60 90 120 150 180 Frequency (Hz) Figure 64: FFT of fluid force response for circular orbit, forward excitation. 0.48 preswirl 84 5.0 Fx 4.0 Fy 3.0 2.0 Reaction Force (N) 1.0 0.0 -1.0 -2.0 -3.0 -4.0 -5.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 65: Monitored force response for multi-frequency elliptical orbit, x excitation. 0.48 preswirl. 1.4 Fx Fy 1.2 1 Reaction Force (N) 0.8 0.6 0.4 0.2 0 0 30 60 90 120 150 180 Frequency (Hz) Figure 66: FFT of fluid force response for elliptical orbit, x excitation. 0.48 preswirl 85 Fx 5.0 Fy 4.0 3.0 2.0 Reaction Force (N) 1.0 0.0 -1.0 -2.0 -3.0 -4.0 -5.0 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) Figure 67: Monitored force response for multi-frequency elliptical orbit, y excitation. 0.48 preswirl. 1.4 Fx 1.2 Fy 1 Reaction Force (N) 0.8 0.6 0.4 0.2 0 0 50 100 150 200 Frequency (Hz) Figure 68: FFT of fluid force response for elliptical orbit, y excitation. 0.48 preswirl 86 5.2.1 Computational time The disadvantage of running a full 360 CFD model is the computational time and resources necessary to achieve a solution. This is particularly the case when using the transient CFD analysis. To realize a practical solution time, parallel processing is necessary and the simulations in this study make use of high-performance computational resources to achieve a faster solution turn around. The simulations were performed using clusters maintained by the institute for cyber- enabled research at Michigan State University. Each cluster consists of similar processors with hardware information detailed in [58]. The steady state eccentric mesh has 8.17 × 106nodes. The meshing strategy was kept the same for the transient mesh, however, with the slightly different concentric geometry, the final node count was 8.72 × 106 nodes. Table 8 shows the simulation times needed to achieve 3 whirl periods for the KCM transient model and one converged steady state solution utilizing 260 cores. The deforming mesh KCM model requires 49.3 hours of simulation time for 4 whirling frequencies and one static eccentric case. Contrarily, the steady state eccentric model requires 5 x 3.5 hrs. = 17.5 hours of computational time. The contrast in computational time between steady state and transient model shows one benefit in using the steady state method when possible. When comparing the KCM and KC transient models, the KCM model provides some computational benefit, as only the forward whirling orbits need to be considered, to achieve an adequate quadratic curve fit, when assuming a circular whirling orbit. Alternatively, the frequency dependent KC model requires 2 linearly independent whirling rotor motions (i.e., forward, and backward), which requires two simulations. In this study, 3 cycles of the base frequency required 60 hours for a combination of 6 excitation frequencies. Therefore, the multifrequency model needed 120 hours to complete 2 simulations. Although time consuming, this method provides less computational effort than running 2 simulations for each whirl frequency 87 individually. Moreover, a multifrequency approach is beneficial when a wide range of excitation frequencies are necessary. Table 8: Transient KCM CFD model computational time Whirl ratio 0% 25% 50% 75% 100% Timesteps/rev -- 360 180 120 90 Periodicity (s) -- 0.024 0.012 0.008 30% Total timesteps -- 1080 540 360 270 Sim. time (hr) 5 21.2 10.65 7.1 5.31 Total (hr) 49.3 5.3. Comparison to Experiment This section compares the methods used to predict the rotordynamic coefficients of the 14 TOR labyrinth seal described in [22]. For the steady state and transient CFD frequency- independent model comparison, the direct and cross-coupled stiffness coefficients from the experiment are determined via the zero-frequency intercept of a second order curve fit of 𝑅𝑒(𝐻𝑖𝑖 ) and 𝑅𝑒(𝐻𝑖𝑗 ), respectively. Moreover, the direct and cross coupled damping are taken from the slope of a linear trendline passing through the origin and fitted to 𝐼𝑚(𝐻𝑖𝑖 ) and 𝐼𝑚(𝐻𝑖𝑗 ), respectively [59]. The frequency dependent KC model is compared to the dynamic coefficients for the entire range of excitation frequency provided in the raw data. Additionally, the frequency dependent KC model considers two rotor whirling models: circular and elliptical. Overall, each CFD methods give a favorable comparison to the experimental rotordynamic terms. In terms of accuracy, the rotordynamic coefficients for each whirling rotor model nearly coincide, despite having different whirling orbits and amplitudes. In other words, these parameters have significant influence on the rotor whirling motions and fluid response forces but have very little influence on the computed rotordynamic coefficients. This is because the peak whirling amplitude is small enough to capture the small perturbation theory represented by the linear reaction force/motion 88 model. Therefore, either model can be used to predict reasonably accurate frequency-dependent rotordynamic coefficients. The rotordynamic coefficients are quantified in Figures 70-85. The largest discrepancy comes when predicting the direct stiffness term. For each swirl case, the frequency independent direct stiffness in Figure 70 is underpredicted and even shows the sign to be negative. This is contrary to experiment which records this value as positive. Arthur [22] reports the same issue when comparing 1 control volume bulk flow model with experimental direct stiffness at medium swirl. However, this term has little influence on whirling rotor stability and although large direct stiffness values can affect the natural frequency and critical speeds, the values recorded by Arthur [22], are too small to influence the critical speeds. Additionally, both the steady state and transient CFD model match the trend with experiment that shows direct stiffness slightly decreasing as pre swirl increases. The multi-frequency transient KC model presents the frequency dependent rotordynamic coefficients. Between CFD predictions and measurements, the frequency-dependent KC CFD model provides a closer agreement to the magnitude of experimental direct damping, particularly at higher frequencies, although the KC model also produces a negative value. Additionally, the multi-frequency KC model predicts the direct stiffness being frequency dependent, especially at higher frequencies, whereas [22] reports the coefficients as frequency independent. However, as noted in [22], the experimental impedances become more erratic at higher frequencies. This partially explains the poor agreement with the frequency independent steady state and transient models. However, it is uncertain why the direct stiffness sign is contrary to experiment. It could be that the direction of the radial forces is defined differently in experiment. 89 Figure 74 shows the cross coupled stiffness term as a function of preswirl ratio. The discrepancy between the steady state and transient CFD predictions is miniscule for all cases, differing by less than 3%. The trend for the CFD shows a linear relationship where the cross coupled stiffness increases as inlet swirl increases which matches the trend in experiment, however, the linear trend shown by the experimental values is more pronounced with a steeper incline. This results in KCM models overpredicting the cross coupled stiffness term at low inlet swirl and under predicting this term at high inlet swirl. The prediction at the medium inlet swirl condition for both steady state and transient CFD gives an almost exact match. Additionally, each model predicts within 25% at the high swirl case. Figure 75-77 shows the frequency dependent cross-coupled stiffness terms for each. In this model, a good match is seen for the cross-coupled stiffness at each pre-swirl ratio over a range of excitation frequencies. However, for the high preswirl case a larger discrepancy is noticed as the CFD underpredicts this value. As mentioned in chapter 4, the rotordynamic coefficients are also influenced by the chosen turbulence model. At this point, further investigation is needed to determine the superiority of turbulence model choice in the multi-frequency whirling rotor transient approach. Nevertheless, the 𝜅-𝜀 model commonly used in industry, is sufficient to provide reasonably accurate cross coupled stiffness predictions with the multifrequency approach. The direct damping term as a function of inlet swirl is shown in Figure 78 for the frequency independent steady state and transient models. The CFD models show a slight increase in damping as inlet swirl increases to the medium inlet swirl case where it begins to decrease as swirl increases. Contrarily, the experimental value shows a continual increase in damping as preswirl increases, however, this increase becomes less pronounced at higher inlet swirl conditions. Although both frequency independent CFD methods underpredict the direct damping, they both come within 90 about 40% of the experimental value. The transient KCM model provides a slightly better direct damping prediction than the steady state model. However, the steady state model maintains preference can be computed in far less time and using less resources. The frequency dependent direct damping terms shown in Figures 79-81 follow the same trend as experiment and provide a close agreement for multiple frequencies. It is also observed that the KC CFD model produces frequency independent results. Lastly, the cross coupled damping coefficient as a function pre swirl is shown in Figure 82. Both CFD models underpredict the cross coupled damping. However, steady state and KCM model predicts a similar trend of increasing cross-coupled damping as preswirl ratio increases. Similar to the direct damping, the cross coupled damping is weakly influenced by the inlet swirl ratios. Moreover, the frequency dependent cross-coupled damping terms shown in Figure 83-85 provide a close agreement with experiment. It is also observed that the KC CFD model produces frequency independent results. 91 2.50 Steady State CFD Transient CFD 2.00 Experiment 1.50 Kxx (MN/m) 1.00 0.50 0.00 -0.50 -1.00 -0.01 0.09 0.19 0.29 0.39 0.49 0.59 Preswirl ratio Figure 69: Comparison of direct stiffness for steady state and transient KCM with experiment. No Swirl CFD, Circle Orbit CFD, Elliptical Orbit 3.0 Experiment 2.5 2.0 1.5 K (MN/m) 1.0 0.5 0.0 -0.5 0 30 60 90 120 150 180 210 -1.0 -1.5 -2.0 Frequency (Hz) Figure 70: Comparison of predicted direct stiffness using two multifrequency whirling models with experiment. 0.08 preswirl. 92 2.5 Medium Swirl 2.0 1.5 1.0 0.5 K (MN/m) 0.0 -0.5 0 30 60 90 120 150 180 210 240 -1.0 -1.5 -2.0 Frequency (Hz) CFD, Circle Orbit CFD, Elliptical Orbit Experiment Figure 71: Comparison of predicted direct stiffness using two multifrequency whirling models with experiment. 0.25 preswirl. 2.0 High Swirl 1.5 1.0 0.5 K (MN/m) 0.0 0 30 60 90 120 150 180 210 -0.5 -1.0 CFD, Circle Orbit -1.5 CFD, Elliptical Orbit Experiment -2.0 Frequency (Hz) Figure 72: Comparison of predicted direct stiffness using two multifrequency whirling models with experiment. 0.48 preswirl. 93 1.40 1.20 1.00 Steady State CFD 0.80 Transient CFD k (MN/m) 0.60 Experiment 0.40 0.20 0.00 0 0.08 0.16 0.24 0.32 0.4 0.48 0.56 Preswirl ratio Figure 73: Comparison of cross coupled stiffness for steady state and transient KCM model with experiment Low Swirl 3.5 CFD, Circle Orbit CFD Elliptical Orbit 2.5 Experiment k (MN/m) 1.5 0.5 -0.5 0 30 60 90 120 150 180 210 -1.5 Frequency (Hz) Figure 74: Comparison of predicted cross coupled stiffness using two multifrequency whirling models with experiment. 0.08 preswirl. 94 Medium Swirl 3.8 3.0 CFD, Circle Orbit CFD, Elliptical Orbit 2.3 Experiment k (MN/m) 1.5 0.8 0.0 0 30 60 90 120 150 180 210 -0.8 -1.5 Frequency (Hz) ` Figure 75: Comparison of predicted cross coupled stiffness using two multifrequency whirling models with experiment. 0.25 preswirl. 2.5 High Swirl 2.0 1.5 k (MN/m) 1.0 0.5 CFD, Circle Orbit CFD, Elliptical Orbit 0.0 Experiment 0 30 60 90 120 150 180 210 Frequency (Hz) Figure 76: Comparison of predicted cross coupled stiffness using two multifrequency whirling models with experiment. 0.48 preswirl. 95 3.00 Steady State CFD Tansient CFD 2.50 Experiment Cxx (kN-s/m) 2.00 1.50 1.00 0.50 -0.01 0.09 0.19 0.29 0.39 0.49 0.59 Preswirl ratio Figure 77: Comparison of direct damping for steady state and transient KCM model with experiment 3.0 Low Swirl 2.6 CFD, Circle Orbit CFD, Elliptical Orbit 2.2 Experiment 1.8 C (kN.s/m) 1.4 1.0 0.6 0.2 -0.2 0 30 60 90 120 150 180 210 Frequency (Hz) Figure 78: Comparison of predicted direct damping using two multifrequency whirling models with experiment. 0.08 preswirl. 96 Medium Swirl 2.7 CFD, Circle Orbit 2.3 CFD, Elliptical Orbit Experiment 1.9 1.5 C (kN.s/m) 1.1 0.7 0.3 -0.1 0 30 60 90 120 150 180 210 -0.5 Frequency (Hz) Figure 79: Comparison of predicted direct damping using two multifrequency whirling models with experiment. 0.25 preswirl. High Swirl 2.7 CFD, Circle Orbit 2.3 CFD, Elliptical Orbit Experiment 1.9 1.5 C (kN.s/m) 1.1 0.7 0.3 -0.1 0 30 60 90 120 150 180 210 240 -0.5 Hz Figure 80: Comparison of predicted direct damping using two multifrequency whirling models with experiment. 0.48 preswirl. 97 Steady State CFD Transient CFD 3.00 Experiment 2.50 cxy (kN-s/m) 2.00 1.50 1.00 0.50 -0.01 0.09 0.19 0.29 0.39 0.49 0.59 Preswirl ratio Figure 81: Comparison of cross coupled damping for steady state and transient KCM model with experiment Low Swirl 2.7 CFD, Circle Orbit 2.3 CFD, Elliptical Orbit Experiment 1.9 1.5 c (kN.s/m) 1.1 0.7 0.3 -0.1 0 30 60 90 120 150 180 210 -0.5 Frequency (Hz) Figure 82: Comparison of predicted cross coupled damping using two multifrequency whirling models with experiment. 0.08 preswirl. 98 Medium Swirl 2.7 CFD, Circle Orbit 2.3 CFD, Elliptical Orbit Experiment 1.9 1.5 c (kN.s/m) 1.1 0.7 0.3 -0.1 0 30 60 90 120 150 180 210 -0.5 Frequency (Hz) Figure 83: Comparison of predicted cross coupled damping using two multifrequency whirling models with experiment. 0.25 preswirl. High Swirl 2.7 CFD, Circle Orbit 2.3 CFD, Elliptical Orbit Experiment 1.9 1.5 c (kN.s/m) 1.1 0.7 0.3 -0.1 0 30 60 90 120 150 180 210 -0.5 Frequency (Hz) Figure 84: Comparison of predicted cross coupled damping using two multifrequency whirling models with experiment. 0.48 preswirl. 99 5.3.1 Stability Characteristics To better understand the effects of the rotordynamic coefficients on the stability of the labyrinth seal, the effective damping is plotted with respect to the excitation frequency. The cross- over frequency at which effective damping is zero indicates stability margin. At excitation frequencies below crossover frequency, the seal is destabilizing. Therefore, it is desirable to have high effective damping ratios. For the 14 TOR labyrinth seal, the CFD in Figures 86-88 shows the crossover frequency increasing as the inlet swirl increases. This is because the destabilizing cross coupled stiffness increases as swirl increases, whereas the stabilizing direct damping is only marginally influenced by the inlet swirl. The trend is the same in the experimental data, however, the crossover frequency occurs at slightly lower frequencies. For low swirl (0.08), the CFD crossover frequency occurs at about 90 Hz. For medium swirl (0.25), the crossover frequency occurs at about 120 Hz. For the high swirl case, the CFD effective damping doesn’t quite enter the stable margin within this range of excitation frequencies. From the experimental data, the effective damping crosses over the stable margin at about 100 Hz and then crosses back over and becomes destabilizing at higher excitation frequencies. This study demonstrates the importance of implementing anti-swirling technology which can help reduce the cross coupling term and improve the range of stable operability. 100 4.0 Low Swirl 3.0 CFD, Circle Orbit CFD, Elliptical Orbit 2.0 Experiment 1.0 Ceff (kN.s/m) 0.0 0 30 60 90 120 150 180 210 -1.0 -2.0 -3.0 -4.0 Frequency (Hz) Figure 85: Comparison of predicted effective damping using two multifrequency whirling models with experiment. 0.08 preswirl. 4.0 Medium Swirl 3.0 CFD, Circle Orbit CFD, Elliptical Orbit 2.0 Experiment 1.0 Ceff (kN.s/m) 0.0 -1.0 0 30 60 90 120 150 180 210 -2.0 -3.0 -4.0 -5.0 Frequency (Hz) Figure 86: Comparison of predicted effective damping using two multifrequency whirling models with experiment. 0.25 preswirl. 101 4.0 High Swirl 3.0 2.0 1.0 0.0 -1.0 0 30 60 90 120 150 180 210 Ceff (kN.s/m) -2.0 -3.0 -4.0 -5.0 CFD, Circle Orbit -6.0 CFD, Elliptical Orbit -7.0 Experiment -8.0 Frequency (Hz) Figure 87: Comparison of predicted effective damping using two multifrequency whirling models with experiment. 0.48 preswirl. 5.4. Conclusions This worked utilized experimental data to compare a steady state and transient CFD approach for predicting rotordynamic coefficients in a long 14-TOR gas labyrinth seal. The first transient approach is a frequency independent KCM model, and the second transient approach transfers the monitored force and displacement data from the time domain to the frequency domain to solve for the frequency dependent stiffness and damping terms. Both the steady state and the frequency independent KCM transient model showed comparable performance and good accuracy predicting the rotordynamic coefficients when compared to experiment. Because the rotordynamic coefficients are about the same between these two methods, the steady state eccentric model is better suited for this analysis, due to its faster computational speed. However, when non- axisymmetric features are present, the KCM model can be applied when the steady state model is no longer applicable. This is the main benefit of using the KCM model, which finds its practicality 102 in seals with non-axisymmetric features, such as: hole pattern seals, seals with anti-swirl vanes, or seals whirling about an eccentric equilibrium [60]. Additionally, the frequency independent KCM model requires less computational time than the multifrequency whirling rotor approach for annular gas seals possessing the same direct force coefficients (𝐾𝑥𝑥 = 𝐾𝑦𝑦 , 𝐶𝑥𝑥 = 𝐶𝑦𝑦 and opposite sign but equal magnitude cross coupling force coefficients (𝑘𝑥𝑦 = −𝑘𝑦𝑥 , 𝑐𝑥𝑦 = −𝑐𝑦𝑥 ). However, the multifrequency whirling rotor models provide superior performance in accuracy and detail in predicting rotordynamic coefficients over a range of excitation frequencies. One benefit is its ability to predict the cross over frequency associated with the effective damping. Nevertheless, due to its long computational times, this model is best suited for labyrinth seals with high running speeds that may exhibit frequency dependent characteristics, particularly at high excitation frequencies. 103 Chapter 6: Transient CFD: Further Investigation 6.1. 18-Tooth on Rotor Balance Piston Seal The 18-tooth on rotor balance piston labyrinth seal operates with an abradable babbitt material, which allows for tighter clearances, while at the same time, mitigates rotordynamic instabilities due to rotor rubs. As the labyrinth teeth rub into the abradable material, a groove is produced, creating a larger area for the fluid to travel through. The groove created can vary in axial width and radial depth. Therefore, the groove depth and width that result from this operation are modeled in this CFD analysis, as opposed to traditional CFD models with a smooth stator. Figure 89 and 90 provide a detailed view of how the balanced piston labyrinth seal is modeled. The abbreviations described in Figures 89 and 90 are defined in Table 9. However, the dimensions are withheld for proprietary reasons. Nevertheless, the ratio of the tooth clearance to the tooth radius 𝐶𝑟 𝐿 is scaled such that = 8.75 × 10−4 and the seal length to diameter ratio is = 0.356. The 𝑟 𝐷 abradable groove parameters are given in Table 10. In this model, it has been observed that under certain conditions, the vortex created in the cavity can change directions. This exhibits one benefit of CFD in capturing the flow characteristics in labyrinth seals that cannot be determined by experiment – namely, which way the flow is spinning. It is predicted that vortex direction will have a noticeable impact on the seal performance. Therefore, it is worthwhile to explore the consequences of this phenomenon, particularly its effect on the rotordynamic coefficients. To accomplish this, this study considers a fixed operating clearance at variable abradable groove widths and depths to cover a scope of groove dimensions. These dimensions are displayed as a ratio of the tooth width (TW), ranging from 20-40% of the tooth width. However, the groove dimensions may not be the only influence on the flow direction 104 within the tooth cavity. The flow rate through the labyrinth seal must also be considered. In this part of the analysis, a constant operating clearance is defined, and the nature of the vortex flow is observed as the compressor goes from choke to surge along the compressor map. Here one can see how/if the flow direction and consequently, the rotordynamic coefficients, are influenced by the maximum flow condition (choke) or minimum flow condition (surge). Figure 88: Balance Piston Labyrinth Seal Dimensions Figure 89: Detailed Labyrinth Seal Tooth Dimensions Table 9: Balance Piston Dimensions TH Tooth height TW Tooth width P Pitch Cr Radial clearance GD Groove dimensions RD Rub dimension 𝛼 Tooth angle 105 Table 10: Abradable groove parametric study Parametric Fixed dimension Variable dimension Axial width RD 30% GD 20% 30% 40% TW TW Radial depth GD 30% RD 20% 30% 40% TW TW The design iterations for a balance piston labyrinth seal are performed using a transient multifrequency approach. The 𝑘−𝜔 based Shear Stress Transport (SST) turbulence model was used to account for fluid turbulence. Furthermore, the CFD simulations are run using fluid properties of a natural gas real gas mixture defined by CoolProp [61].The balance piston seal is defined as one fluid domain, which also includes P2 inject sitting above the 5th tooth. The seal domain is given a specified rotating speed and the boundaries are defined using pressures and temperatures as shown in Figure 13. The static pressure from the back face of the impeller cavity provide the inlet pressure conditions into the balance piston. A proper grid density for the 18-tooth on rotor balance piston labyrinth seal is also analyzed in this section. Unlike the other seals, this seal is not experimental, but is used in the industrial centrifugal compressors of Solar Turbines Inc. Consequently, this labyrinth seal is largest seal and operates at the highest pressure conditions. Therefore, the mesh size can easily become excessive if not properly optimized. Hence, the principals established in chapter 3.5 are implemented for the 18-tooth on rotor labyrinth seal. A section cut view of the balance piston seal is shown in Figure 91. Additionally, the most important influence on the mesh is the circumferential node. Figure 92 shows the circumferential node count mesh independents for the 𝐾𝑥𝑥 and 𝑘𝑥𝑦 for the following design point. 228 circumferential nodes provides the optimal mesh size for accuracy and computational time. 106 Figure 90: Balance piston seal mesh, section cut view 35000 45000 kxy 40000 Kxx 30000 35000 25000 30000 kxy (lbf/in) Kxx (lbf/in) 20000 25000 15000 20000 15000 10000 10000 5000 5000 0 0 100 125 150 175 200 225 250 275 300 325 Circumferential nodes Figure 91: Circumferential node study for 18 TOR balance piston seal 6.2. Flow Characteristics This section explores the flow characteristics of the balance piston labyrinth seal at various operating conditions and tooth clearances. Figure 93a-c shows pressure distribution across the labyrinth seal with P2 inject for three flow points: choke, design point, and surge. This figure shows that the pressure distribution across the labyrinth seal is similar for each flow condition. Because the static pressure increases through the return vane, the pressure entering from P2 inject is the highest at this point, forcing a portion of the flow back toward the rear cavity of the last 107 impeller stage. This amounts to roughly 1/3 of the leakage through the labyrinth seal travelling back toward the rear cavity. The remaining 2/3 of the leakage exits the balance piston outlet and is recirculated back toward suctions. As mentioned in section 2.1, flow will travel through tooth clearance, expand in the cavity, and takes on a circumferential nature. However, there is no way to test which way the flow is circulating within the cavity, i.e., clockwise (CW, or counterclockwise (CCW). (a) (b) (c) Figure 92: Pressure distribution through the labyrinth seal. a) choke b) midpoint c) surge Therefore, CFD becomes a powerful tool in understanding the flow pattern through the labyrinth seal. For a smooth TOS or TOR seal (no abradable grooves), the jet flow will typically attach to 108 Figure 93: Detailed example of traditional flow pattern through labyrinth seal the wall parallel to the tooth clearance exit and create a vortex when a portion of the flow attaches to the cavity wall prior to exiting the cavity. This flow pattern was portrayed Figure 3 of section 2.2. and is reproduced in more detail in Figure 94. Moreover, this flow pattern is commonly represented in the literature [37, 41, 49]. However, with abradable grooves and P2 inject in the TOR labyrinth seal, there is added complexity to the flow pattern. This research shows that under certain conditions, the vortex created in the tooth cavity can change directions (e.g., clockwise to counterclockwise or vice versa), as the flow promulgates through the seal. To demonstrate this, a velocity vector profile at 15300 rpm is shown in Figure 95 for a full 360-degree transient model with circular whirl. One can notice that the flow changes directions from clockwise to counterclockwise shortly after P2 inject. This occurs at the 4th and 9th cavity away from P2 inject, moving towards the labyrinth seal outlet back toward suction. This flow behavior has only scarcely been reported in the literature [62], and to this authors knowledge, has not been investigated for a balance piston labyrinth seal with P2 inject. Moreover, this phenomenon is not detected when running the steady state model. This demonstrates a major advantage of using transient CFD over 109 steady-state CFD techniques which can provide a more detailed fluid flow field beneficial in understanding the development of fluid excitation forces in gas labyrinth seals. Figure 94: Velocity vectors in stationary frame Taking a closer look, the abradable grooves can create four areas of recirculation as depicted in Figure 96. The first area of recirculation occurs above the leading tooth edge where part of the jet stream flow travels through the tooth clearance where it runs into a portion of the flow deflected upward from the tooth tip (zone 1). Similarly, a vortex is created in the upper right corner of the abradable groove before exiting the tooth (zone 2). A small recirculation zone also occurs right above the tooth tip (zone 3). Lastly, when the tooth is close enough to the babbitt groove, a nozzle is formed by the groove vertical edge and the side surface of the tooth. As the flow accelerates through the nozzle, the flow becomes detached from the shroud wall. Consequently, an area of recirculation occurs in this region. This creates another vortex in the cavity, labeled zone 4. For comparison, the same operating conditions of the design point are implemented on a theoretical balance piston labyrinth seal without grooves, with the tooth clearance remaining the same. Figure 97 shows the flow pattern in this configuration behaves similar to a conventional seal portrayed in Figure 94, i.e., the jet flow attaches to the stationary wall, and is carried over to the next tooth, creating a clockwise vortex in each cavity. Moreover, in this scenario, the four areas of recirculation are not realized. This small but important detail has 110 a major influence on the flow characteristics and leakage of the seal. The circumferential velocity profile in the stationary frame shown in Figure 98 gradually gains speed as it passes through each seal cavities, creating a positive swirl component. In contrast, the area increase between the rotating and stationary components with the Babbitt groove creates a more turbulent flow path which acts as a de-swirling mechanism, reducing the circumferential velocity component. This contrast is portrayed in the stationary frame circumferential velocity contour plot in Figure 98 and quantified in the swirl coefficient in Figure 100. Consequently, as the circumferential velocity decreases, the meridional velocity component slightly increases. This, along with the area increase contributes to the mass flow rate increasing in the grooved seal. The leakage comparison for each scenario is displayed in Table 11. The % indicates the fraction of the total leakage that exits the labyrinth seal outlet and the fraction that is directed back toward the rear cavity of the last compressor stage. In both scenarios, about 2/3rds of the leakage exits the labyrinth seal outlet and 1/3rd travels back toward the rear cavity. Moreover, there is only a marginal difference in the pressure drop through each cavity caused by the abradable grooves. Since the circumferential velocity component has a major impact on the stability characteristics of the seal, quantified in the cross coupled stiffness term. This could provide some benefit in reducing the destabilizing cross coupled stiffness term. However, the recirculation zones defined in Figure 96 and the presence of flow reversal caused by the abradable grooves, introduce additional mechanisms that may influence the stability caused by the fluid excitation forces in labyrinth seals. 111 Figure 95: Recirculation zones for TOR seal with abradable grooves Figure 96: Flow pattern for TOR seal with smooth stator Figure 97: Circumferential velocity in the stationary frame, no abradable grooves 112 Figure 98: Circumferential velocity in the stationary frame, with abradable grooves Rear P2 inject Seal outlet 0.80 Cavity 0.70 0.60 0.50 Swirl 0.40 0.30 0.20 With grooves Without grooves 0.10 0.00 0 0.2 0.4 0.6 0.8 1 Axial Location (x/L) Figure 99: Swirl coefficient, with and without abradable grooves Figure 100: Meridional velocity in the stationary frame, no abradable grooves 113 Figure 101: Meridional velocity in the stationary frame, abradable grooves Table 11: Leakage comparison, groove, with and without abradable grooves Grooves No Grooves 𝑙𝑏 𝑙𝑏 Leakage ( ) % Leakage ( ) % 𝑠 𝑠 Inject 4.82 3.21 Outlet -3.10 64.3% -2.14 66.6% Inlet (opening) -1.72 35.6% -1.12 33.4% Additionally, Figures 103-114 display the balance piston labyrinth seal velocity profiles as the compressor operates from choke to surge. The nature of the flow pattern is marginally affected by the operating conditions changing from choke to surge. The quantified swirl coefficient in Figure 116 is about the same at each flow condition. Moreover, the flow direction behaves the same for each flow condition, changing from clockwise to counterclockwise in the 4th and 9th cavity. Figure 102: Circumferential velocity in the stationary frame, surge 114 Figure 103: Circumferential velocity in the stationary frame, midpoint (c) Figure 104: Circumferential velocity in the stationary frame, choke Figure 105: Meridional velocity, surge Figure 106: Meridional velocity, midpoint 115 Figure 107: Meridional velocity, choke Figure 108: Velocity vectors in stationary frame, surge, 4th cavity CCW flow 116 Figure 109: Velocity vectors in stationary frame, surge, 9th cavity CCW flow Figure 110: Velocity vectors in stationary frame, midpoint, 4th cavity CCW flow Figure 111: Velocity vectors in stationary frame, midpoint, 9th cavity CCW flow 117 Figure 112:Velocity vectors in stationary frame, choke, 4th cavity CCW flow Figure 113: Velocity vectors in stationary frame, choke, 9th cavity CCW flow 118 Rear Seal outlet P2 inject 0.80 Cavity 0.70 0.60 0.50 Swirl 0.40 0.30 0.20 Surge Midpoint 0.10 Choke 0.00 0 0.2 0.4 0.6 0.8 1 Axial Location (x/L) Figure 114: Swirl coefficient, choke to surge 6.4. Rotordynamic Coefficients This section focuses on identifying the rotordynamic coefficients for the 18 tooth on rotor balance piston seal using the transient multifrequency CFD approach with a circular orbit. The excitation frequency range is given in Table 12, which corresponds to a whirling frequency of 0.30 ×, 0.60 ×, 0.90 ×, 1.20 ×, and 1.50 ×, the running speed. Table 12: Rotor whirling frequency Frequency (Hz) 77 153 230 306 383 Whirl speed (rad/s) 481 961 1442 1923 2403 % Running speed 30% 60% 90% 120% 150% First, a comparison is made between a case with abradable grooves to one without abradable grooves. Additionally, rotordynamic coefficients along three points on a compressor map are analyzed and compared, demonstrating how the coefficients change along the compressor characteristics. First, Figures 115-119 show frequency dependent rotordynamic coefficients for the C31 balance piston seal with and without abradable grooves. Case 1 represents the balance 119 piston seal with modeled grooves and case 2 is the model with a smooth stator. By comparison, the impedance forces are larger when the grooves aren’t present. This was predicted when it was discovered that the circumferential velocity component was higher when abradable grooves are not modeled However, this discrepancy between case 1 and case 2 become more pronounced at higher frequencies. Figure 115 demonstrates the frequency dependent direct stiffness terms for each case. At subsynchronous excitation frequencies (< 1 × 𝑟𝑢𝑛𝑛𝑖𝑛𝑔 𝑠𝑝𝑒𝑒𝑑), there is marginal difference in the direct stiffness between each case. As the excitation frequency increases beyond synchronous whirl (1 × 𝑟𝑢𝑛𝑛𝑖𝑛𝑔 𝑠𝑝𝑒𝑒𝑑), the magnitude of the direct stiffness value also increases. This increase is greater when the abradable grooves are not present. Therefore, the abradable grooves display a large impact on the frequency dependency of the seal. Interaction between the acoustic waves in a circumferential annulus and the rotor excitation frequency is the cause of the predicted frequency dependency [63]. The cross coupled stiffness values for each case are represented in Figure 116. As predicted, the cross-coupled stiffness values are higher when considering the case with a smooth stator. Similar to the direct stiffness, the difference between the two cases is smaller at low excitation frequencies. However, this difference becomes more pronounced after passing through an excitation frequency of around 0.9 × 𝑡ℎ𝑒 𝑟𝑢𝑛𝑛𝑖𝑛𝑔 𝑠𝑝𝑒𝑒𝑑, or 230 Hz. At this point, the cross coupled stiffness of the labyrinth seal with the smooth stator increases at a higher rate and becomes more destabilizing. This shows that the abradable grooves influence the frequency dependency of the seal by reducing it. Figure 116 shows the frequency dependent direct damping terms which increase in a linear fashion as the excitation frequency increases. Although the abradable grooves supply some rotordynamic benefit by reducing the cross-coupled stiffness, this benefit is slightly counteracted by a reduction in the stabilizing direct damping term. The smooth stator in case 2 provides more damping than the seal with abradable 120 grooves. This relationship is quantified by considering the effective damping for both cases displayed in Figure 118. The higher damping available from the geometry without grooves is 𝑟𝑎𝑑 enough to maintain similar effective damping up until a whirl frequency of 1442 𝑠 (0.9𝑥 𝑟𝑢𝑛𝑛𝑖𝑛𝑔 𝑠𝑝𝑒𝑒𝑑) or 230 Hz. Moreover, the effective damping is positive up to a 𝑟𝑎𝑑 subsynchronous whirl frequency of 961.32 (0.6 × 𝑟𝑢𝑛𝑛𝑖𝑛𝑔 𝑠𝑝𝑒𝑒𝑑) or 153 Hz, which is 𝑠 stabilizing. The effective damping then crosses over to a negative value, contributing to a destabilizing characteristic. This decline is more dramatic for case 2 (without grooves) as the whirl frequency approaches synchronous excitation and especially in the super synchronous range (> 1 × 𝑟𝑢𝑛𝑛𝑖𝑛𝑔 𝑠𝑝𝑒𝑒𝑑). Lastly, the cross-coupled damping terms in Figure 119 are higher when the abradable grooves are not modeled. One can also see an increasing relationship between the cross-coupled stiffness and the excitation frequency. This is not the case when abradable grooves are present, which again demonstrates the influence of the abradable grooves have on reducing the frequency dependency of the rotordynamic coefficients. 121 0 -50,000 0 500 1000 1500 2000 2500 3000 -100,000 K (Case 1) -150,000 K (Case 2) -200,000 K (lbf/in) -250,000 -300,000 -350,000 -400,000 -450,000 -500,000 Whirl Frequency (rad/s) Figure 115: Direct stiffness coefficient 550,000 475,000 400,000 k (Case 1) k (Case 2) 325,000 k (lbf/in) 250,000 175,000 100,000 25,000 -50,000 0.0 500.0 1000.0 1500.0 2000.0 2500.0 3000.0 Whirl Frequency (rad/s) Figure 116: Cross coupled stiffness coefficient 122 140 120 100 C (Case 1) C (lbf-s/in) 80 C (Case 2) 60 40 20 0 0 500 1000 1500 2000 2500 3000 Whirl Frequency (rad/s) Figure 117: Direct Damping 150 -50 0 500 1000 1500 2000 2500 3000 -250 Ceff (lbf-s/in) -450 -650 Ceff (Case 1) Ceff (Case 2) -850 -1,050 -1,250 Whirl Frequency (rad/s) Figure 118: Effective damping 123 120 100 80 c (Case 1) c (lbf-s/in) 60 c (Case 2) 40 20 0 0 500 1000 1500 2000 2500 3000 Whirl Frequency (rad/s) Figure 119: Cross-coupled damping Figures 120-124 shows the frequency dependent rotordynamic coefficients for the balance piston labyrinth seal at different operating conditions between choke and surge with abradable grooves modeled. In this case, speed remains the same, so the main influencing parameter is the pressure ratio and flow rate. First, the direct stiffness is show in Figure 120 shows that the magnitude of the direct stiffness increases from choke to surge. Moreover, the direct stiffness shows frequency dependency at lower frequency and is highly frequency dependent after 150 Hz (0.6 × 𝑟𝑢𝑛𝑛𝑖𝑛𝑔 𝑠𝑝𝑒𝑒𝑑). Although cross coupled stiffness is usually considered frequency independent, this Figure 121 shows the k values to be slightly frequency dependent at lower frequencies and are highly frequency dependent after 223 Hz (0.9𝑥 𝑟𝑢𝑛𝑛𝑖𝑛𝑔 𝑠𝑝𝑒𝑒𝑑). Additionally, at subsynchronous whirling frequencies, the destabilizing cross coupled stiffness is highest at choked flow. As the fluid conditions approach synchronous whirl and beyond, the cross coupled stiffness is higher at conditions near surge. After 223 Hz, the cross coupled stiffness at 124 each flow point are found to be positive and increase with increasing frequency. Nevertheless, this destabilizing term is counterbalanced by an increase in direct damping as the compressor goes from choke to surge. Moreover, this term also demonstrates frequency dependence as it increases linearly with excitation frequency. The effective damping stability characteristic is shown in Figure 123. This figure demonstrates that there is no considerable difference in the labyrinth seal stability based on the operating flow conditions. As the compressor goes along the compressor map, the operating conditions only play a minor role in influencing the stability of the machine. Lastly, the cross coupled damping predicted in Figure 124 shows the cross coupled damping increasing from choke to surge. 0 -30,000 K (Surge) -60,000 K (Midpt) -90,000 K (Choke) -120,000 K (lbf/in) -150,000 -180,000 -210,000 -240,000 -270,000 -300,000 0 500 1000 1500 2000 2500 3000 Whirl Frequency (rad/s) Figure 120: Direct stiffness from choke to surge 125 210,000 180,000 150,000 k (Surge) 120,000 k (Midpt) k (lbf/in) 90,000 k (Choke) 60,000 30,000 0 -30,000 0 500 1000 1500 2000 2500 3000 -60,000 Whirl Frequency (rad/s) Figure 121: Cross coupled stiffness from choke to surge 60 50 40 C (lbf-s/in) C (Surge) 30 C (Midpt) C (Choke) 20 10 0 0 500 1000 1500 2000 2500 3000 Whirl Frequency (rad/s) Figure 122: Direct damping from choke to surge 126 175 Ceff (Surge) Ceff (Midpt) 100 Ceff (Choke) 25 -50 0 500 1000 1500 2000 2500 3000 Ceff (lbf-s/in) -125 -200 -275 -350 -425 -500 Whirl Frequency (rad/s) Figure 123: Effective damping from choke to surge 35 30 25 c (lbf-s/in) 20 c (Surge) c (Midpt) 15 c (Choke) 10 5 0 0 500 1000 1500 2000 2500 3000 Whirl Frequency (rad/s) Figure 124: Cross coupled damping from choke to surge 127 6.4.1 Abradable grooves Up to this point, it has been demonstrated that the abradable grooves associated with the 18 tooth on rotor labyrinth seal have a major impact on the leakage flow pattern through the seal. Additionally, the geometrical parameters of the abradable grooves become an additional influencer to the rotordynamic behavior of the machine. Therefore, this study concludes by examining the effects of the abradable groove parameters possible during operation of a centrifugal compressor. First, the rub dimension is fixed (see Figure 90), and the groove width dimension is altered to capture the effects of the groove width on the rotordynamic coefficients. In the next series of cases, the tooth groove width is fixed, and the rub dimension is altered for three cases to understand the effects of rub depth on the rotordynamic coefficients. First, the flow characteristics are compared for the balance piston seal with various axial groove widths. For the abradable groove case with an axial width of 40% of the TW, the counterclockwise flow direction in the 4th and 9th cavity after P2 inject present in the 30% axial groove dimension reverts back to a clockwise rotation. This flow pattern remains clockwise through the rest of the seal. Reducing the axial groove width to 20% creates a scenario where the flow direction 50% of the cavities to the right of P2 inject, experience counterclockwise rotation and 50% undergo clockwise rotation. This influence on the swirl velocity is demonstrated in Figure 126-127 which shows the swirl ratio through the balance piston labyrinth seal as the axial width dimension changes. Additionally, the 3 dimensional swirl pattern in Figures 128-129 demonstrates the contrast between the axial groove widths of 20% and 40%. Finally, the axial groove width also affects the areas of recirculation as defined in Figure 93. Figures 130-131 demonstrate this comparison between the axial groove width of 20% and 40%. As the axial groove width increases, the recirculation defined as zone 4 in Figure 96 becomes smaller. 128 Figure 125: Circumferential velocity in the stationary frame, 20% axial groove width Figure 126: Circumferential velocity in the stationary frame, 40% axial groove width Figure 127: 3D streamlines in tooth cavity, 20% in. axial width 129 Figure 128: 3D streamlines in tooth cavity, 40%. axial width Rear P2 inject Seal outlet 0.80 Cavity 0.70 0.60 0.50 Swirl 0.40 0.30 0.20 Axial width 40% Axial width 30% 0.10 Axial width 20% 0.00 0 0.2 0.4 0.6 0.8 1 Axial Location (x/L) Figure 129: Swirl coefficient, axial groove width comparison 130 Figure 130: Labyrinth tooth recirculation zone with reduced axial width, 20% in. Figure 131: Labyrinth tooth recirculation zone with increased axial width, 40%. Figures 133-137 shows the influence of the axial groove width on the rotordynamic coefficients. In Figure 133, the direct stiffness remains negative for each case and increases in magnitude as the axial width increases from 20% to 40% on each side of the labyrinth tooth. In the subsynchronous excitation range, the cross coupled stiffness is highest at the smallest axial groove width. This is somewhat counterintuitive because the swirl component significantly reduces when the flow 131 direction changes from clockwise to counterclockwise. Therefore, the change in flow direction to counterclockwise seems to promote an increase in the destabilizing cross coupled stiffness term. As the excitation frequency passes synchronous, the widest axial groove width contributes to a higher cross coupled stiffness term. Moreover, the rotordynamic comparison of the labyrinth seal with and without grooves demonstrated that the groove dimensions have an influence on the frequency dependency of the seal, particularly at higher frequencies. The destabilizing 𝑘 factor at low excitation frequencies is slightly counteracted by an increase in the direct damping (Figure 132) as the axial groove width increases, producing a similar effective damping at each excitation frequency. At higher frequencies, the larger axial groove width becomes slightly more destabilizing and the difference between each case becomes more noticeable. However, at this point, each case is operating below the stable margin. Lastly, the cross coupled damping also increases as the axial width increases. 0 -50,000 K (Axial width 40%) K (Axial width 30%) -100,000 K (Axial width 20%) K (lbf/in) -150,000 -200,000 -250,000 -300,000 0 500 1000 1500 2000 2500 3000 Whirl Frequency (rad/s) Figure 132: Direct stiffness, axial groove width comparison 132 210,000 180,000 150,000 120,000 k (lbf/in) 90,000 k (Axial width 40%) 60,000 k (Axial width 30%) 30,000 k (Axial width 20%) 0 -30,000 -60,000 0 500 1000 1500 2000 2500 3000 Whirl Frequency (rad/s) Figure 133: Cross coupled stiffness, axial groove width comparison 60 50 40 C (lbf-s/in) 30 C (Axial width 40%) 20 C (Axial width 30%) C (Axial width 20%) 10 0 0 500 1000 1500 2000 2500 3000 Whirl Frequency (rad/s) Figure 134: Direct damping, axial groove width comparison 133 200 Ceff (Axial width 40%) Ceff (Axial width 30%) 100 Ceff (Axial width 20%) 0 0 500 1000 1500 2000 2500 3000 Ceff (lbf-s/in) -100 -200 -300 -400 -500 Whirl Frequency (rad/s) Figure 135: Effective damping, axial groove width comparison 35 30 25 c (lbf-s/in) 20 15 c (Axial width 40%) c (Axial width 30%) 10 c (Axial width 20%) 5 0 0 500 1000 1500 2000 2500 3000 Whirl Frequency (rad/s) Figure 136: Cross coupled damping, axial groove width comparison 134 Next, the flow characteristics are compared for the balance piston seal with various radial groove depths. A parametric study of the effects of the radial groove depth shows that this parameter also has an important influence on the flow pattern within the seal. This influence on the swirl velocity is demonstrated in Figure 138-140 which shows the swirl ratio through the balance piston labyrinth seal as the radial width dimension changes. At a radial groove depth of 20% the TW, the flow pattern rotates clockwise in each cavity to the right of P2 inject, which is similar to the case examining an axial groove width of 40%. As the radial depth dimension decreases, the recirculation zones initially defined in Figure 96 are reduced in size, and the fluid leaving the nozzle formed by is formed by the groove vertical edge and the side surface of the tooth exits in a more axial direction, allowing the jet stream flow to stay attached to the stator surface. Figure 137: Circumferential velocity in the stationary frame, 20% in. radial depth Figure 138: Circumferential velocity in the stationary frame, 40%. radial depth 135 Rear P2 inject Seal outlet 0.80 Cavity 0.70 0.60 0.50 Swirl 0.40 0.30 0.20 Radial depth 40% Radial depth 30% 0.10 Radial depth 20% 0.00 0 0.2 0.4 0.6 0.8 1 Axial Location (x/L) Figure 139: Swirl coefficient, radial groove width comparison Figure 140: Labyrinth tooth recirculation zone with reduced radial depth, 20% in. 136 Figure 141: Labyrinth tooth recirculation zone with increased radial depth, 40%. Figure 143 displays the direct stiffness as this parameter increases from 20% to 40%. There is no clear trend relating the abradable groove depth with the direct stiffness. At subsynchronous whirl frequencies, the largest radial groove depth (40%) predicts the highest magnitude of direct stiffness. The lowest radial groove depth model predicts slightly higher magnitudes of direct stiffness than the medium groove depth of 30%. In Figure 144, the radial groove depth seems to influence the cross coupled stiffness by increasing it as the radial depth increases. At this parameter, the fluid path experiences the most flow reversals from clockwise to counterclockwise rotation in the cavity. Similar flow characteristics were present when examining the labyrinth seal with a minimum abradable groove width of 20%. This case also had a higher a destabilizing cross coupled stiffness term, particularly at subsynchronous excitation frequencies. This indicates that the flow direction plays an important role in influencing the seals stability characteristics. The direct damping is only marginally effected by the radial groove depth. Therefore, the smallest 137 radial groove depth (20%) predicts a slightly higher effective damping. Lastly, Figure 147 shows that the cross coupled damping generally increases as the radial depth decreases. Additionally, the rub groove depth dimension has a slight influence on the frequency dependence of the cross coupled damping term, as this coefficient linearly increases much faster as the excitation frequency increases. 0 -50,000 K (Radial depth 40%) K (Radial depth 30%) -100,000 K (Radial depth 20%) K (lbf/in) -150,000 -200,000 -250,000 -300,000 0 500 1000 1500 2000 2500 3000 Whirl Frequency (rad/s) Figure 142: Direct stiffness, radial depth comparison 138 210,000 180,000 150,000 120,000 k (lbf/in) 90,000 k (Radial depth 40%) 60,000 k (Radial depth 30%) 30,000 k (Radial depth 20%) 0 -30,000 0 500 1000 1500 2000 2500 3000 -60,000 Whirl Frequency (rad/s) Figure 143: Cross coupled stiffness, radial depth comparison 60 50 40 C (lbf-s/in) 30 20 C (Radial depth 40%) C (Radial depth 30%) 10 C (Radial depth 20%) 0 0 500 1000 1500 2000 2500 3000 Whirl Frequency (rad/s) Figure 144: Direct damping, radial depth comparison 139 200 Ceff (Radial depth 40%) Ceff (Radial depth 30%) 100 Ceff (Radial depth 20%) 0 Ceff (lbf-s/in) -100 -200 -300 -400 -500 0 500 1000 1500 2000 2500 3000 Hz Figure 145: Effective damping, radial depth comparison 35 30 25 c (lbf-s/in) 20 15 c (Radial depth 40%) c (Radial depth 30%) 10 c (Radial depth 20%) 5 0 0 500 1000 1500 2000 2500 3000 Whirl Frequency (rad/s) Figure 146: Cross couped damping, radial depth comparison 140 6.5. Conclusions Experimental test seals reported in the literature are characterized as seal-only models and represent scaled representatives of most industrial seals. This work provides the benefit of analyzing a full scale balance piston labyrinth seal used in industry. A transient CFD model was used to explore the flow characteristics of the seal and to determine the frequency dependent rotordynamic coefficients. Moreover, the CFD was modeled with and without abradable grooves, and conclusions can be drawn between the two cases. The abradable grooves cause an increase in the leakage, which increases as the rub radial depth, or the axial groove width increases. When comparing the rotordynamic terms, the seal with abradable grooves predicted slightly lower cross coupled stiffness and direct damping terms. Therefore, the effective damping was about the same for each case at subsynchronous excitation frequencies. However, the seal with a smooth stator has a higher frequency dependency and increases more rapidly at higher frequencies. It was also discovered that the geometrical parameters associated with the abradable grooves have a major influence on the flow pattern within the seal that can’t be overlooked. If the labyrinth teeth rub deeper into the abradable material on the stator, the resultant flow path can cause the flow to change directions (CW to CCW) within the cavity. The same is true when the axial groove width is kept to a minimum dimension. This flow direction has an important impact on the rotordynamic terms. When the labyrinth seal experience more CCW flow within the cavities, the cross coupled stiffness term tended to be higher. Lastly, the rotordynamic coefficients were predicted at three flow points along the compressor map, to determine how these stability characteristics behave from choke to surge. There is only a slight variation in the flow pattern through the labyrinth seal at each condition and the difference in rotordynamic coefficients for each case is marginal, particularly at low excitation frequencies. 141 142 Chapter 7: Conclusions 7.1. Summary and Recommendations Computational fluid dynamics has made great strides over the years in prediction fluid- dynamic and rotordynamic performance in annular gas labyrinth seals used in centrifugal gas compressors. This is demonstrated in the CFD work presented here. First, the steady state CFD approach was examined. This approach gives favorable results when applied to labyrinth seals with axisymmetric geometries. Moreover, when compared to the transient CFD approach, the steady state model offers great superiority in computational speed and therefore is still relevant for rotordynamic analysis today. However, when seal geometries become non-axisymmetric, a transient model is necessary. In this case, a KCM model which assumes frequency independent coefficients can be used with good accuracy. Moreover, this approach provides reasonable computational times depending on the number of discrete whirling frequencies necessary to achieve an appropriate curve fit. Both the steady state and the frequency independent KCM transient model showed comparable performance in predicting the rotordynamic coefficients when compared to experiment. Lastly, a multifrequency whirling rotor approach was also considered. The difference between this model and the transient KCM model is that the predicted rotordynamic coefficients are frequency dependent. Recent studies have shown that under certain conditions, labyrinth seals can possess frequency dependent rotordynamic coefficients. In the multifrequency orbit approach, a linear superposition of several whirling frequency excitations is used to determine the rotordynamic coefficients over a range of frequencies. A prescribed circular orbit or elliptical orbit can be used, and the influence of orbital shape (circle or elliptical) is insignificant in estimating the rotordynamic coefficients from CFD. The transient CFD approach applied to the balance piston labyrinth seal offered valuable insight in the flow mechanisms that influence the rotordynamic coefficients. It was found that the 143 grooves created in the abradable material have a significant impact on the fluid flow characteristics in the seal. When the radial rub depth is large, or the axial groove width is small, the fluid direction within the labyrinth seal cavities can change direction from CW to CCW. The discovery of this phenomenon gives greater insight on the influences that the flow characteristics have on rotordynamic performance. When the labyrinth seal experience more CCW flow within the cavities, the cross coupled stiffness term tended to be higher. This work only briefly touched on the influence turbulence model using the steady state approach. It may be beneficial to expand upon this work and explore the superiority of turbulence models using a transient approach for labyrinth seal applications. Additionally, the balance piston labyrinth seal in this work only considered one speed line. Further investigation is needed to determine the influence of flow conditions over a scope of rotational speeds. Additionally, although the transient CFD approach provides great benefit in accuracy and detail, the scope of its application still remains hindered by extra computational time and effort. Therefore, before expanding the scope of this study to other seal geometries and operating conditions, it would be beneficial to reduce the computational time for the transient solver. A few possibilities to do this are: to decrease the mesh size or increase the time step. Of these, the mesh size in this study was established from a grid independent study which determined the optimum mesh size for speed and accuracy. However, each mesh was created in ANSYS Mesh and therefore it is possible that other meshing software’s may provide superior performance in meshing quality and speed, which is worth looking into. Increasing the number of time steps without sacrificing accuracy will need further investigation 144 REFERENCES 145 REFERENCES [1] D. G. Wilson, “The design of high-efficiency turbomachinery and gas turbines.”, MIT Press, 2014. [2] H. Krain, "Review of Centrifugal Compressor’s Application and Development," Journal of Turbomachinery-transactions of The Asme, 2005. [3] K. H. Lüdtke, Process Centrifugal Compressors, Springer-Verlag Berlin Heidelberg, 2004. [4] A. Engeda, "From the Crystal Palace to the Pump Room," Mechanical Engineering, pp. 50- 53, 1999. [5] H. Bloch, Compressors And Modern Process Applications, JOHN WILEY & SONS, INC, 2001. [6] H. Bloch, A practical guide to compressor technology, JOHN WILEY & SONS, INC, 1995. [7] J. Stahley, Dry Gas Seals Handbook, Tulsa, OK: PennWell Publishing, 2005. [8] T. Gresh, Compressor Performance: Aerodynamics for the User, Elsevier Science & Technology Books, 2001. [9] A. R. H., "Mean Streamline Aerodynamic Performance Analysis of Centrifugal Compressors," Journal of Turbomachinery , vol. 117, pp. pp. 360-366, 1995. [10] B. Qiao, Y. Ju and C. Zhang, "Numerical Investigation on Labyrinth Seal Leakage Flow and Its Effects on Aerodynamic Performance for a Multistage Centrifugal Compressor," Journal in Fluids Engineering of the ASME, vol. 141, 2019. [11] R. Marechale, M. Ji and M. Cave, "EXPERIMENTAL AND NUMERICAL INVESTIGATION OF LABYRINTH SEAL CLEARANCE IMPACT ON CENTRIFUGAL IMPELLER PERFORMANCE," in Proceedings of ASME Turbo Expo 2015: Turbine Technical Conference and Exposition, Montréal, Canada, 2015. [12] R. Flitney, "Chapter 3: Rotary Seals," in Seals and Sealing Handbook, Sixth Edition) ed., ScienceDirect, 2014, pp. 105-288. [13] S. Yoon and P. Allaire , Control of Surge in Centrifugal Compressors by Active Magnetic Bearings, Springer, 2012. [14] N. Rieger, Rotordynamics 2 - Problems in Turbomachinery, Springer-Verlag, 1988. 146 [15] H. B. Jeffcott, "The Lateral Vibration of Loaded Shafts in the Neighborhood of a Whirling Speed The Effect of Want of Balance," Phil. Lofagazine Series 6, vol. 37, p. 304, 1919. [16] J. M. Vance, Machinery Vibration and Rotordynamics, New York: Wiley, 2010. [17] K. Kwanka, "Dynamic Coefficients of Stepped Labyrinth Gas Seals," Journal of Engineering for Gas Turbines and Power, pp. pp. 473-477, 2000. [18] D. Childs, Turbomachinery Rotordynamics: Phenomena, Modeling, and Analysis, New York: John Wiley & Sons, 1993. [19] A. Picardo and D. Childs, "Rotordynamic Coefficients for a Tooth-on-Stator Labyrinth Seal at 70 Bar Supply Pressures: Measurements Versus Theory and Comparisons to a Hole Pattern Seal," Journal of Engineering for Gas Turbines and Power, vol. 127, pp. 843-855, 2005. [20] P. Hanlon, Compressor Handbook, New York: McGraw-Hill Companies, Inc., 2001. [21] A. A. Fox, "An Examination of Gas Compressor Stability and Rotating Stall," Solar Turbines Inc. , San Diego, CA. [22] P. Arthur, "Measured Rotordynamic Coefficients and Leakage for a Tooth-on-Rotor Labyrinth Seal with Comparisons to a Tooth-on-Stator Labyrinth Seal," M.S.M.E. Thesis, Department of Mechanical Engineering, Texas A&M University., 2015. [23] A. Picardo, "High Pressure Testing of See-Through Labyrinth Seals," M.S. Thesis, Department of Mechanical Engineering, Texas A&M University, 2003. [24] N. Mehta, "Comparison of a Slanted-Tooth See-Through Labyrinth Seal to a Straight- Tooth See-Through Labyrinth Seal for Rotordynamic Coefficients and Leakage," MS Thesis, Department of Mechanical Engineering, Texas A&M University, College Station, 2012. [25] J. S. Alford, "Protecting Turbomachinery from Self-Excited Rotor Whirl," ASME. J. Eng. Power, 1965. [26] Wachter Joachim and H. Benckert, "Flow Induced Spring Coefficients of Labyrinth Seals for Application in Rotor Dynamics," NASA Conf. Pub. 2133, Vols. Rotordynamic Instability Problems in High-Performance Turbomachinery, Texas A&M Workshop, pp. 189-212, 1980. [27] R. Kirk, "Evaluation of Aerodynamic Instability Mechanisms for Centrifugal Compressors – Part 2," Journal of Vibrations, Acoustic, Stress, and Reliability in Design, vol. 110, pp. 207-212, 1988. 147 [28] J. Pelletti, "A Comparison of Experimental Results and Theoretical Predictions for the Rotordynamic Coefficients of Short (L/D=1/6) Labyrinth Seals," M.S.M.E. Thesis, Texas A&M University , 1990. [29] J. K. Scharrer, "A Comparison of Experimental and Theoretical Results for Labyrinth Gas Seals," Doctoral Dissertation, Texas A & M University, College Station, TX, 1987. [30] T. Iwatsubo, "Evaluation Of Instability Forces Of Labyrinth Seals In Turbines Or Compressors," Proc. Rotordynamic Instability Problems In High Performance Turbomachinery, Nasa, pp. 139-167, 1980. [31] H. R. Wyssmann, T. C. Pham and R. J. Jenny, "Prediction Of Stiffness And Damping Coefficients For Centrifugal Labyrinth Seal," Journal Of Engineering For Gas Turbines And Power, pp. Pp. 920-926, 1984. [32] D. Childs and J. Scharrer, "An Iwatsubo-Based Solution for Labyrinth Seals: Comparison to Experimental Results," Journal of Engineering for Gas Turbines and Power, vol. 108, pp. 599-604, 1986. [33] B. Williams and R. Flack, "Calulations of Rotor Dynamic Coefficients for Labyrinth Seals," International Journal of Rotating Machinery, vol. Vol. 4, no. No. 4, pp. 257-269, 1998. [34] D. Rhode , S. Hensel and M. Guidry, "Labyrinth Seal Rotordynamic Forces Using 3- Dimenstional Navier Stokes Code," ASME J. Tribology, pp. 683-689, 1992. [35] J. Moore, "Three-dimensional CFD rotordynamic analysis of gas labyrinth seals," Journal of Vibrations and Acoustics, vol. 125, pp. 427-433, 2003. [36] M. Athavale , R. C. Hendricks and A. J. Przekwas, "A 3D-CFD code for accurate prediction of fluid flows and fluid forces in seals.," in Proceedings of the Advanced ETO Propulsion Conference, NASA CP3282, NASA Marshall Space Flight Center, Huntsville, AL, 1994. [37] T. Hirano, Z. Guo and R. G. Kirk, "Application of computational fluid dynamics analysis for rotating machinery—part II- labyrinth seal analysis," in Transactions of the ASME, 2005. [38] R. G. Kirk and Z. Guo, "Calibration of Labyrinth Seal Bulk Flow Design Analysis Predictions to CFD Simulation Results," in Proc. IMECHE 8th International Conference on Vibrations in Rotating Machinery, 2004. [39] E. D. Thompson, "Study of Forces and Dynamic Coefficients in Whirling and Whirling and Eccentric Labyrinth Seals," M.S. Thesis Virginia Tech University , 2009. 148 [40] M. Wensheng , L. Zhusheng and C. Zhao-bo, "Study of leakage and stability in labyrinth seal using CFD," in Proceedings of 2011 International Conference on Electronic & Mechanical Engineering and Information Technology 588-591, 2011. [41] R. Gao, "Computational Fluid Dynamic and Rotordynamic Study on the Labyrinth Seal," PhD Dissertation Virginia Polytechnic Institute and State University, 2012. [42] S. Subramanian, A. Sekhar and B. Prasad, "Rotordynamic characteristics of rotating labyrinth gas turbine seal with centrifugal growth," Tribology International 97, pp. 349- 359, 2016. [43] T. Wu and L. San Andres, "Leakage and Dynamic Force Coefficients for Two Labyrinth Gas Seals: Teeth-on-Stator and Interlocking Teeth Configurations. A Computational Fluid Dynamics Approach to Their Performance," ASME. J. Eng. Gas Tur, 2018. [44] G. Chochua and T. Soulas, "Numerical Modeling of Rotordynamic Coefficients for Deliberately Roughened Stator Gas Annular Seals," Journal of Tribology-transactions of The Asme, pp. 424-429, 2006. [45] X. Yan and Z. Feng, "Investigations on the rotordynamic characteristics of a hole-pattern seal using transient CFD and periodic circular orbit model," ASME Journal of Vibration and Acoustics, vol. 133 (4), 2011. [46] Z. Li, J. Li and Z. Feng, "Numerical comparisons of rotordynamic characteristics for three types of labyrinth gas seals with inlet preswirl," Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 2016. [47] L. Zhigang, J. Li and Z. Feng, "Comparisons of Rotordynamic Characteristics Predictions for Annular Gas Seals Using the Transient Computational Fluid Dynamic Method Based on Different Single-Frequency and Multifrequency Rotor Whirling Models," Journal of tribology - Transaction of ASME, 2016. [48] K. Kwanka, "Rotordynamic Coefficients of Short Labyrinth Gas Seals—General Dependency on Geometric and Physical Parameters," Tribology Transactions, pp. 558-563, 2007. [49] G. Wagner , K. Steff, R. Gausmann and M. Schmidt, "Investigation on the Dynamica Coefficients of Impeller Eye Labyrinth Seal," in Proceedings of the 38th Turbomachinery Symposium, 2009. [50] U. Bolleter, A. Wyss and R. Sturchler, "Measurement of Hydrodynamic Interaction Matrices of Boiler Feed Pump Impellers," ASME Journal of Vibrations, Stress and Reliability in Design, vol. 109, pp. 144-151, 1987. 149 [51] H. Ulbrich, "New Techniques Using Magnetic Bearings," in Magnetic Bearings, New York, Springer Verlag, 1988, pp. 281-288. [52] R. Gao and G. Kirk, "CFD Study on Stepped and Drum Balance Labyrinth Seal," Tribology Transactions, vol. 56, pp. 663-671, 2013. [53] "Turbulence and Wall Function Theory," in ANSYS CFX-Theory Solver Guide, Canonsburg, PA, ANSYS INC., 2011, pp. pp. 89-154. [54] P. A. Davidson, "Turbulence: An Introduction for Scientists and Engineers," in Turbulent shear flows and simple closure models, Oxford, Oxford University Press, 2015, pp. 107- 191. [55] D. Wilcox, Turbulence Modeling for CFD, La Cãnada, CA: DCW Industries, Inc,, 1993. [56] F. Menter, "Review of the shear-stress transport turbulence model experience from an industrial perspective," International Journal of Computational Fluid Dynamics, vol. 23, pp. 305-316, 2009. [57] M. R. Thorat and J. R. Hardin, "Rotordynamic Characteristics Prediction for Hole-Pattern Seals Using Computational Fluid Dynamics," ASME. J. Eng. Gas Turbines Power, vol. 142, 2020. [58] A. Keen and C.-M. Chang, "2018 Cluster Resources," Michigan State University, 2020. [Online]. Available: https://wiki.hpcc.msu.edu/pages/viewpage.action?pageId=20120131. [Accessed 2021]. [59] S. Aurthur and D. Childs, "MEASURED ROTORDYNAMIC AND LEAKAGE CHARACTERISTICS OF A TOOTH-ON-ROTOR LABYRINTH SEAL WITH COMPARISONS TO A TOOTH-ON-STATOR LABYRINTH SEAL AND Predictions," in Proceedings of ASME Turbo Expo 2015: Turbine Technical Conference and Exposition, Montreal, 2015. [60] K. Hoopes, J. Moore, A. Rimpel, C. Kulhanek and B. Venkataraman, "A METHOD FOR ROTORDYNAMIC FORCE PREDICTION OF A CENTRIFUGAL COMPRESSOR IMPELLER FRONT CAVITY USING A TRANSIENT WHIRLING CFD TECHNIQUE," in Proceedings of ASME Turbo Expo 2019, Phoenix, Arizona, 2019. [61] I. H. Bell, J. Wronski, and S. Quoilin, "Pure and Pseudo-pure Fluid Thermophysical Property Evaluation and the Open-Source Thermophysical Property Library CoolProp," Industrial & Engineering Chemistry Research, vol. 53, no. 6, 2014. 150 [62] D. Zhang, C. Lee and M. Cave, "A CFD STUDY ON THE DYNAMIC COEFFICIENTS OF LABYRINTH SEALS," in Proceedings of ASME Turbo Expo 2012, Copenhagen , 2012. [63] V. Giuseppe, M. R. Thorat, D. W. Childs and M. Libraschi, "IMPACT OF FREQUENCY DEPENDENCE OF GAS LABYRINTH SEAL ROTORDYNAMIC COEFFICIENTS ON CENTRIFUGAL COMPRESSOR STABILITY," in Proceedings of ASME Turbo Expo 2010: Power for Land, Sea and Air, Glasgow, UK, 2010. [64] J. K. Whalen, E. (. Alvarez and L. P. Palliser, "Thermoplastic Labyrinth Seals for Centrifugal Compressors," in Whalen, John K.; Alvarez, Eduardo (Ed); Palliser, Lester P. (2004). Thermoplastic Labyrinth Seals for Centrifugal Compressors. Texas A&M University. Turbomachinery Laboratories, Texas A&M, 2004. [65] R. Tiwari, Theory & Practice of Rotor Dynamics, Indian Institute of Technology Guwahati. [66] U. Yèucel, "Leakage and swirl velocities in labyrinth seals," M.S. Thesis Lehigh University, 1996. [67] "Turbulence and Wall Function Theory," in ANSYS CFX-Solver Theory Guide, Canonsburg, PA, ANSYS inc., 2011, pp. ch. 2, pp. 89-154.. 151