UFDC Home  myUFDC Home  Help 



Full Text  
COMPUTATIONAL MODELING OF THERMODYNAMIC EFFECTS IN CRYOGENIC CAVITATION By YOGEN UTTURKAR A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2005 Copyright 2005 by Yogen Utturkar To my wife and parents ACKNOWLEDGMENTS I am grateful to several individuals for their support in my dissertation work. The greater part of this work was made possible by the instruction of my teachers, and the love and support of my family and friends. It is with my heartfelt gratitude that I acknowledge each of them. Firstly, I would like to express sincere thanks and appreciation to my advisor, Dr. Wei Shyy, for his excellent guidance, support, trust, and patience throughout my doctoral studies. I am very grateful for his remarkable wisdom, thoughtprovoking ideas, and critical questions during the course of my research work. I thank him for encouraging, motivating, and always prodding me to perform beyond my own limits. Secondly, I would like to express sincere gratitude towards my coadvisor, Dr. Nagaraj Arakere, for his firm support and caring attitude during some difficult times in my graduate studies. I also would like to express my appreciation to the members of my dissertation committee Dr. Louis Cattafesta, Dr. James Klausner, and Dr. Don Slinn, for their valuable comments and expertise to better my work. I deeply thank Dr. Siddharth Thakur (ST) for providing me substantial assistance with the STREAM code and for his cordial suggestions on research work and career planning. My thanks go to all the members of our lab, with whom I have had the privilege to work. Due to the presence of all these wonderful people, work is more enjoyable. In particular, it was a delightful experience to collaborate with Jiongyang Wu, Tushar Goel, and Baoning Zhang on various research topics. I would like to express my deepest gratitude towards my family members. My parents have always provided me unconditional love. They have always given top priority to my education, which made it possible for me to pursue graduate studies in the United States. I would like to thank my grandparents for their selfless affection and loving attitude during my early years. My wife's parents and her sister's family have been extremely supportive throughout my graduate education. I greatly appreciate their trust in my abilities. Last but never least, I am thankful beyond words to my wife, Neeti Pathare. Together, we have walked through this memorable, cherishable, and joyful journey of graduate education. Her honest and unfaltering love has been my most precious possession all these times. I thank her for standing besides me every time and every where. To Neeti and my parents, I dedicate this thesis! TABLE OF CONTENTS page A C K N O W L E D G M E N T S ................................................................................................. iv LIST OF TABLES ....................................................... ............ ....... ....... ix LIST OF FIGURES ............................... .... ...... ... ................. .x LIST OF SY M B OL S .... .................................................. .. ....... ............. xiv A B S T R A C T ...............x.................v.................................................. xv iii CHAPTER 1 INTRODUCTION AND RESEARCH SCOPE ........................................................1 1.1 T ypes of C aviation ................ ...2..... ... ................ ........................ ...2 1.2 Cavitation in Cryogenic Fluids Thermal Effect................... ...............4 1.3 Contributions of the Current Study.................... ...................8 2 LITER A TU RE R EV IEW ............................................................. ....................... 9 2.1 General Review of Recent Studies ........................................ ...... ............... 9 2.2 M odeling Therm al Effects of Cavitation ......................................... .................19 2.2.1 Scaling Law s ................................. ... .. ..... ........... ... 19 2.2.2 E xperim ental Studies.............................................. ............... ... 24 2.2.3 Numerical Modeling of Thermal Effects ................................................27 3 STEADY STATE COM PUTATION S............................ ....................................... 33 3.1 G governing E quations .................................................. .............................. 33 3.1.1. C aviation M modeling ..................................................... ...................35 3.1.1.1 M erkle et al. M odel ................................. ...... .............. 35 3.1.1.2 Sharp Interfacial Dynamics Model (IDM) ..................................... 35 3.1.1.3 Mushy Interfacial Dynamics Model (IDM) ...................................37 3.1.2 Turbulence M odeling ............................................................ ...............40 3.1.3 Speed of Sound (SoS) M odeling .... .......... ...................................... 41 3.1.4 T herm al M odeling ......................................................................... ... ... 42 3.1.4.1 Fluid property update ......................... ............ ... .................42 3.1.4.2 Evaporative cooling effects............... .............................................43 3.1.5 B oundary C onditions...................................................................... .. .... 44 3.2 R results and D discussion ................................................. .............................. 45 3.2.1 Cavitation in N oncryogenic Fluids ................................. ............... 45 3.2.2 Cavitation in Cryogenic Fluids........................................ ............... 50 3.2.2.1 Sensitivity analyses ........................ ....... ..............52 3.2.2.2 Assessment of cryogenic cavitation models over a wide range of c o n d itio n s ...................................... ...................................... 5 8 4 TIMEDEPENDENT COMPUTATIONS FOR FLOWS INVOLVING PHASE CH AN GE ............................................................... ..... ..... ......... 69 4 .1 G alliu m F u sion ................................................................... 7 1 4 .1.1 G governing E quations ...................................................................... .. .... 72 4.1.2 N um erical A lgorithm ...................................................................... .. .... 73 4 .1.3 R esu lts ...................... ....... ......................... ................. 7 8 4.1.3.1 A accuracy and grid dependence ................................ ............... 79 4 .1.3 .2 S tab ility ............... .... ..... ... .. .............................. 82 4.1.3.3 Data analysis by reducedorder description ...................................83 4.2 Turbulent Cavitating Flow under Cryogenic Conditions ...................................88 4.2.1 G governing Equations ........................................................ ............... 88 4.2.1.1 Speed of sound modeling .............. ............................................ 89 4.2.1.2 Turbulence m odeling.................................... ........................ 89 4.2.1.3 Interfacial velocity m odel...................................... ............... 89 4.2.1.4 B oundary conditions ............................................. ............... 91 4.2.2 N um erical A lgorithm ........................................... ........................... 91 4 .2 .3 R e su lts ................................................................9 4 5 SUMMARY AND FUTURE WORK ............................................ ............... 101 5 .1 S u m m a ry ............................................................................................1 0 1 5.2 Future W ork ................................. .............................. ........ 103 APPENDIX A BACKGROUND OF GLOBAL SENSITIVITY ANALYSIS..............................105 B REVIEW AND IMPLEMENTATION OF POD ................................................107 B .1 R ev iew ....................................................... ................. 107 B .2 M them atical B background ...................................................... ..... .......... 110 B .3 N um erical Im plem entation ...................................... ......................... ........... 112 B.3.1 Singular Value D ecom position (SVD ) ...................................................112 B .3.2 Postprocessing the SVD Output........................................................... 113 B .4 F low chart .................................................................. .... ......................... 114 B.5 Computing Efforts...... .................. ........ ........ 114 L IST O F R E FE R E N C E S ......................................................................... ................... 116 BIOGRAPHICAL SKETCH ............................................................. ..................126 LIST OF TABLES Table pge 1. Properties of some cryogens in comparison to water at N.B.P, 1.01 bars...................6 2. Source term s in cavitation m models ...................................................... ......... ......... 13 3. V ariants of the k m odel..................... ........................................... ............... 16 4. Summary of studies on thermal effects in cavitation................... .................32 5. Flow cases chosen for the hydrofoil geometry. .................................. .................51 6. Flow cases chosen for the ogive geometry. ...................................... ............... 51 7. Location of the primary vortex for the St = 0.042, Ra = 2.2 x 105 and Pr = 0.0208 c a se ...................................... ..................................................... 8 0 LIST OF FIGURES Figure p 1. Different types of cavitation (a) Traveling cavitation (b) Cloud cavitation (c) Sheet cavitation (d) Supercavitation (e) Vortex cavitation.......................................3 2. Saturation curves for water, Nitrogen, and Hydrogen............................... ...............5 3. Phasic densities along liquidvapor saturation line for water and liquid Nitrogen.........7 4. General classification of numerical methods in cavitation....................................11 5. Variation of Speed of Sound with phase fraction ......... ............................ 14 6. Two cavitation cases for Bfactor analysis ....................................... ............... 21 7. Schematic of bubble model for extracting speed of sound........................ ........27 8. Schematic of cavity models (a) Distinct interface with vaporous cavity (Sharp IDM) (b) Smudged interface with mushy cavity (Mushy IDM).............................36 9. Behavior of p,/ p and p, /p vs. a, for the two models; pl/v = 100 and / = 0 .0 9 .............................................................................4 0 10. Pressuredensity and pressureenthalpy diagrams for liquid nitrogen in the liquid vapor saturation regime (Lemmon et al. 2002). Lines denote isotherms in K elv in ............................................................................... 4 3 11. Illustration of the computational domains for hemispherical projectile and NACA66MOD hydrofoil (noncryogenic cases) ..........................................46 12. Pressure coefficients over the hemispherical body (o = 0.4); D is the diameter of the hemispherical projectile. (a) Impact of grid refinement for Mushy IDM (b) Comparison between pressure coefficients of different models on the coarse grid.47 13. Cavity shapes and flow structure for different cavitation models on hemispherical projectile (o = 0.4). (a) Merkle et al. Model (b) Sharp IDM (c) Mushy IDM........48 14. Pressure coefficient over the NACA66MOD hydrofoil at two different cavitation numbers; D is hydrofoil chord length. (a) a = 0.91(b) a = 0.84 ..........................49 15. Illustration of the computational domain accounting the tunnel for the hydrofoil (Hord 1973a) and 0.357inch ogive geometry (Hord 1973b) (cryogenic cases)......50 16. Noncavitating pressure distribution (a) case '290C', D represents hydrofoil thickness and x represents distance from the circular bend (b) case '312D', D represents ogive diameter and x represents distance from the leading edge...........52 17. Sensitivity of Merkle et al. Model prediction (surface pressure and temperature) to input parameters namely Cdet and Cprod for the hydrofoil geometry....................53 18. Main contribution of each design variable to the sensitivity of Merkle et al. (1998) model prediction; case '290C' (a) Surface pressure (b) Surface temperature .........55 19. Pressure and temperature prediction for Merkle et al. Model for the case with best match with experimental pressure; Ce,, = 0.85; t = 0.85; = 1.1; = 0.9 ...........56 20. Pressure and temperature prediction for Merkle et al. Model for the case with best match with experimental temperature; Ct = 1.15; t ; = 0.85; p =1.1; L = 1.1 ........56 21. Sensitivity of Mushy IDM prediction for case '290C' (surface pressure and temperature) to the exponential transitioning parameter P....................................58 22. Surface pressure and temperature for 2D hydrofoil for all cases involving liquid Nitrogen. The results referenced as 'Mushy IDM' and 'Merkle et al. Model' are contributions of the present study. ........................................ ....................... 59 23. Cavitation number (o = p p py(T)/(0.5pU )) based on the local vapor pressure Merkle et al. Model. Note the values of ao = p p p(TJ)/(0.5pU2 ) for the cases '290C' and '296B' are 1.7 and 1.61, respectively. .............. ................... 61 24. Cavity shape indicated by liquid phase fraction for case '290C'. Arrowed lines denote streamlines (a) Merkle et al. Model isothermal assumption (b) Merkle et al. Model with thermal effects (c) Mushy IDM with thermal effects ............62 25 Evaporation (ti ) and condensation (i+ ) source term contours between the two cavitation models case '290C'. Refer to equations (3.8) and (3.18) for the form ulations. .........................................................................63 26. Surface pressure and temperature for 2D hydrofoil for cases involving liquid Hydrogen. The results referenced as 'Mushy IDM' and 'Merkle et al. Model' are contributions of the present study. ........................................ ....................... 65 27. Surface pressure and temperature for axisymmetric ogive for all the cases (Nitrogen and Hydrogen). The results referenced as 'Mushy IDM' and 'Merkle et al. M odel' are contributions of the present study.................................................67 28. Schematic of the 2D Gallium square geometry with the Boundary Conditions (S h y y et al. 19 9 8 ) ................................................................... 7 1 29. 2D interface location at various instants for St = 0.042, Ra = 2.2 x 105 and Pr = 0.0208. White circles represent interface locations obtained by Shyy et al. (1995) on a 41x41 grid at time instants at t = 56.7s, 141.8s, & 227s respectively. .79 30. Grid sensitivity for the St = 0.042, Ra = 2.2 x 105 and Pr = 0.0208, 2D case (a) Centerline vertical velocity profiles at t = 227s (b) Flow structure in the upper left domain at t = 57s; 41x41 grid (c) Flow structure in the upperleft domain at t = 57s; 81x81 grid. ......................................................................80 31. 3D interface location at t = 57s & 227s for St = 0.042, Ra = 2.2 x 105 and Pr = 0.0208 case on a 41x41x41 grid. Top and bottom: adiabatic; North and West: T = 0; South and East: T = 1 (heated walls) ..................................... ............... 82 32. Interface location and flow pattern for the 3D case, St = 0.042, Ra = 2.2 x 105 and Pr = 0.0208 case at various z locations, at t = 227s...................... ...................82 33. POD modes showing velocity streamlines (y /(r); i = 1, 2,3,4 ) for St = 0.042, Ra = 2.2 x 103 and Pr = 0.0208 case. q(r,t) = (r,t). ......................................... 84 34 Scalar coefficients (q (t);i = 1, 2,..., 8) for St = 0.042, Ra = 2.2 x 103 and Pr = 0 .0208 case. q (r,t) = V (r,t) ..................................... ...................... ............... 85 35. Time instants when coefficients of respective POD modes show the first peak; St = 0.042, Ra = 2.2 x 103 & Pr = 0.0208 case......................................... ..............86 36. Horizontalcenterline vertical velocities for St = 0.042, Ra = 2.2 x 103 and Pr = 0 .02 0 8 case at t = 3 s ................................................................. 86 37. Comparison between CFD solution and reconstructed solution for St = 0.042, Ra = 2.2 x 103 and Pr = 0.0208 case at t = 20s. Contours represent velocity magnitude and lines represent streamlines ........... .............................................87 38. Comparison between CFD solution and reconstructed solution for St = 0.042, Ra = 2.2 x 103 and Pr = 0.0208 case at t = 3s. Contours represent velocity magnitude and lines represent streamlines ........... .............................................88 39. Pressure history at a point near inlet, A = 0.11 D, where D is the chord length of the hydrofoil (a) '283B (b) '290C' .......................................................................... 95 40. Cavity shape and flow structure for case '283B'. Lines denote streamlines .............95 41. Pressure history at a point near the inlet for the case '290C'. St denotes the non dimensional perturbation frequency (fD / U, = ft, )................................ 96 42. Pressure history at a point near the inlet for the cases '296B' and '290C'. The nondimensional perturbation frequency (fD/U, = ft.) is 1.0..............................97 43 POD modes for the velocity field; V, (r); i = 1, 2,3; q(r, t) = V(r, t) ...........................98 44 POD modes for the velocity field; f, (r); i = 1, 2,3; q(r, t) = V(r, t) ...........................98 LIST OF SYMBOLS c: cavitation number 3: boundary layer/cavity thickness L: length (of cavity or cavitating object) q: heat flux; generalized flow variable in POD representation V: total volume t: time x, y, z: coordinate axes r: position vector S: streamwise direction in the curvilinear coordinate system i, k, n: indices u, v, w: velocity components p: pressure T: temperature p: density a: volume fraction f mass fraction h: sensible enthalpy s: entropy c: speed of sound 7: ratio of the two specific heats for gases r: stress tensor P: production of turbulent energy C1, CE2, ok, oE, C, C: constants Q: total kinetic energy (inclusive of turbulent fluctuations) k: turbulent kinetic energy E : turbulent dissipation F: filter function for filterbased modeling / : dynamic viscosity K: thermal conductivity v: volume flow rate L: latent heat C : specific heat C,: pressure coefficient a: thermal diffusivity ih : volume conversion rate B: Bfactor to gauge thermal effect E: dimensional parameter to assess thermal effect ,7: control parameter in the Mushy IDM; coefficient of thermal expansion g: gravitational acceleration U: velocity scale D: length scale (such as hydrofoil chord length or ogive diameter) R: bubble radius CFL: Courant, Freidricks, and Levy Number St: Stefan Number Ra: Rayleigh Number Pr: Prandtl Number Re: Reynolds Number A: difference; filter size in filterbased turbulence modeling V: gradient or divergence operator y/: POD mode 0: timedependent coefficient in the POD series r : energy content in the respective POD mode Subscripts: oo : reference value (typically inlet conditions at the tunnel) 0: initial conditions s: solid 1: liquid v: vapor m: mixture f friction Q: discharge c: cavity L: laminar t: turbulent I: interfacial n: normal to local gradient of phase fraction sat: saturation conditions dest: destruction of the phase prod: production of the phase +: condensation : evaporation A, H, S, M, G, B: terms in a discretized equation nb: neighboring nodes P: at the cell of interest Superscripts/Overhead symbols: +: condensation : evaporation ': fluctuating component *: normalized value; updated value in the context of PISO algorithm : vector : average : Favreaveraged n: time step level k: iteration level Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy COMPUTATIONAL MODELING OF THERMODYNAMIC EFFECTS IN CRYOGENIC CAVITATION By Yogen Utturkar August 2005 Chair: Wei Shyy Cochair: Nagaraj Arakere Major Department: Mechanical and Aerospace Engineering Thermal effects substantially impact the cavitation dynamics of cryogenic fluids. The present study strives towards developing an effective computational strategy to simulate cryogenic cavitation aimed at liquid rocket propulsion applications. We employ previously developed cavitation and compressibility models, and incorporate the thermal effects via solving the enthalpy equation and dynamically updating the fluid physical properties. The physical implications of an existing cavitation model are reexamined from the standpoint of cryogenic fluids, to incorporate a mushy formulation, to better reflect the observed "frosty" appearance within the cavity. Performance of the revised cavitation model is assessed against the existing cavitation models and experimental data, under noncryogenic and cryogenic conditions. Steady state computations are performed over a 2D hydrofoil and an axisymmetric ogive by employing real fluid properties of liquid nitrogen and hydrogen. The thermodynamic effect is demonstrated under consistent conditions via the reduction in xviii the cavity length as the reference temperature tends towards the critical point. Justifiable agreement between the computed surface pressure and temperature, and experimental data is obtained. Specifically, the predictions of both the models are better; for the pressure field than the temperature field, and for liquid nitrogen than liquid hydrogen. Global sensitivity analysis is performed to examine the sensitivity of the computations to changes in model parameters and uncertainties in material properties. The pressurebased operator splitting method, PISO, is adapted towards typical challenges in multiphase computations such as multiple, coupled, and nonlinear equations, and sudden changes in flow variables across phase boundaries. Performance of the multiphase variant of PISO is examined firstly for the problem of gallium fusion. A good balance between accuracy and stability is observed. Timedependent computations for various cases of cryogenic cavitation are further performed with the algorithm. The results show reasonable agreement with the experimental data. Impact of the cryogenic environment and inflow perturbations on the flow structure and instabilities is explained via the simulated flow fields and the reduced order strategy of Proper Orthogonal Decomposition (POD). CHAPTER 1 INTRODUCTION AND RESEARCH SCOPE The phenomenon by which a liquid forms gasfilled or vaporfilled cavities under the effect of tensile stress produced by a pressure drop below its vapor pressure is termed cavitation (Batchelor 1967). Cavitation is rife in fluid machinery such as inducers, pumps, turbines, nozzles, marine propellers, hydrofoils, journal bearings, squeeze film dampers etc. due to wide ranging pressure variations along the flow. This phenomenon is largely undesirable due to its negative effects namely noise, vibration, material erosion etc. Detailed description of these effects can be readily obtained from literature. It is however noteworthy that the cavitation phenomenon is also associated with useful applications. Besides drag reduction efforts (Lecoffre 1999), biomedical applications in drug delivery (Ohl et al. 2003) and shock wave lithotripsy (Tanguay and Colonuis 2003), environmental applications for decomposing organic compounds (Kakegawa and Kawamura 2003) and water disinfection (Kalamuck et al. 2003), and manufacturing and material processing applications (Soyama and Macodiyo 2003) are headed towards receiving an impetus from cavitation. Comprehensive studies on variety of fluids such as water, cryogens, and lubricants have provided significant insights on the dual impact of cavitation. Experimental research methods including some mentioned above have relied on shock waves (Ohl et al. 2003), acoustic waves (Chavanne et al. 2002), and laser pulses (Sato et al. 1996), in addition to hydrodynamic pressure changes, for triggering cavitation. Clearly, research potential in terms of understanding the mechanisms and characteristics of cavitation in different fluids, and their applications and innovation is tremendous. The applicability and contributions of the present study within the above mentioned framework are described later in this chapter. 1.1 Types of Cavitation Different types of cavitation are observed depending on the flow conditions and fluid properties. Each of them has distinct characteristics as compared to others. Five major types of cavitation have been described in literature. They are as follows: (a) Traveling cavitation It is characterized by individual transient cavities or bubbles that form in the liquid, expand or shrink, and move downstream (Knapp et al. 1970). Typically, it is observed on hydrofoils at small angles of attack. The density of nuclei present in the upstream flow highly affects the geometries of the bubbles (Lecoffre 1999). Traveling cavitation is illustrated in Figure 1(a). (b) Cloud cavitation It is produced by vortex shedding in the flow field and is associated with strong vibration, noise, and erosion (Kawanami et al. 1997). A reentrant jet is usually the causative mechanism for this type of cavitation (Figure l(b)). (c) Sheet cavitation It is also known as fixed, attached cavity, or pocket cavitation (Figure l(c)). Sheet cavitation is stable in quasisteady sense (Knapp et al. 1970). Though the liquidvapor interface is dependent on the nature of flow, the closure region is usually characterized by sharp density gradients and bubble clusters (Gopalan and Katz 2000). (d) Supercavitation Supercavitation can be considered as an extremity of sheet cavitation wherein a substantial fraction of the body surface is engulfed by the cavity (Figure l(d)). It is observed in case of supersonic underwater projectiles, and has interesting implications on viscous drag reduction (Kirschner 2001). (e) Vortex cavitation It is observed in the core of vortices in regions of high shear (Figure 1(e)). It mainly occurs on the tips of rotating blades and in the separation zone of bluff bodies (Knapp et al. 1970). Figure 1. Different types of cavitation* (a) Traveling cavitation (b) Cloud cavitation (c) Sheet cavitation (d) Supercavitation (e) Vortex cavitation Reproduced from Franc et al. (1995) with permission from EDP Sciences (b) 1.2 Cavitation in Cryogenic Fluids Thermal Effect Cryogens serve as popular fuels for the commercial launch vehicles while petroleum, hypergolic propellants, and solids are other options. Typically, a combination of liquid oxygen (LOX) and liquid hydrogen (LH2) is used as rocket propellant mixture. The boiling points of LOX and LH2 under standard conditions are 183 F and 423 F, respectively. By cooling and compressing these gases from regular conditions, they are stored into smaller storage tanks. The combustion of LOX and LH2 is clean since it produces water vapor as a byproduct. Furthermore, the power/gallon ratio of LH2 is high as compared to other alternatives. Though storage, safety, and extreme low temperature limits are foremost concerns for any cryogenic application, rewards of mastering the use of cryogens as rocket propellants are substantial (NASA Online Facts 1991). A turbopump is employed to supply the low temperature propellants to the combustion chamber which is under extremely high pressure. An inducer is attached to the turbopump to increase its efficiency. Design of any space vehicle component is always guided by minimum size and weight criteria. Consequently, the size constraint on the turbopump solicits high impellor speeds. Such high speeds likely result in a zone of negative static pressure (pressure drop below vapor pressure) causing the propellant to cavitate around the inducer blades (Tokumasu et al. 2002). In view of the dire consequences, investigation of cavitation characteristics in cryogens, specifically LOX and LH2, is an imperative task. Intuitively, physical and thermal properties of a fluid are expected to significantly affect the nature of cavitation. For example, Helium4 shows anomalous cavitation properties especially past the Xpoint temperature mainly due to its transition to superfluidity (Daney 1988). Besides, quantum tunneling also attributes to cavity formation in Helium4 (Lambare et al. 1998). Cavitation of Helium4 in the presence of a glass plate (heterogenous cavitation) has lately produced some unexpected results (Chavanne et al. 2002), which are in contrast to its regular cavitating pattern observed under homogenous conditions (without a foreign body). Undoubtedly, a multitude of characteristics and research avenues are offered by different types of cryogenic fluids. The focus of the current study is, however, restricted to cryogenic fluids such as LOX, LH2, and liquid Nitrogen due to their aforesaid strong relevance in space applications. It is worthwhile at the outset to contrast their behavior/physical properties to water. Water Nitrogen "50 i Temperture (K) Tempewatre (K) lo Hydrogen I 7 225 5 "250' 27S" 0 5 3 Temperature (K) Figure 2. Saturation curves for water, Nitrogen, and Hydrogen Obtained from REFPROP v 7.0 by Lemmon et al. (2002) Table 1. Properties of some cryogens in comparison to water at N.B.P, 1.01 bars Substance Specific heat Liquid Liquid/Vapor Thermal Vaporization (J/Kg.K) density density conductivity heat (kg/m3) (W/mK) (KJ/Kg) Water 4200 958 1603 681 2257 H2 9816 71 53 100 446 N2 2046 809 175 135 199 02 1699 1141 255 152 213 Source: Weisend et al. (1998) Refer to Figure 2 and Table 1. The operating point of cryogenic liquids is generally quite close to the critical point unlike water. Furthermore, as indicated by Figure 2, the saturation pressure curves for cryogens demonstrate a much steeper slope v/s temperature, as compared to water. Consequently, the vapor pressure of liquids such as LH2 and liquid Nitrogen is expected to show great sensitivity to small temperature drops. Since cavitation is predominantly governed by the vapor pressure, the significance of this thermodynamic sensitivity on the flow problem is clear upfront. Further insight can be obtained from Table 1. Liquidtovapor density ratio in the case of cryogenic fluids is substantially lower than water. Thus, these fluids require a greater amount of liquid, and in turn latent heat, than water to cavitate (form vapor) under similar flow conditions. Furthermore, the thermal conductivity for the low temperature fluids is consistently lower than water. These facts indicate that the sensiblelatent heat conversion in cryogenic fluids is expected to develop a noticeable temperature gradient surrounding the cavitation region. The impact of this local temperature drop is magnified by the steep saturation curves observed in Figure 2. Subsequently, the local vapor pressure experiences a substantial drop in comparison to the freestream vapor pressure leading to suppression in the cavitation intensity. Experimental (Hord 1973a, 1973b) and numerical results (Deshpande et al. 1997) on sheet cavitation have shown a 2040% reduction in cavity length due to the thermodynamic effects in cryogenic fluids. Water Liq. N2 ....... ............ ... ... .. .. TemperaluM wK) g g S~ i E J1 s 0 A i asa ,, 7 J s : Figure 3. Phasic densities along liquidvapor saturation line for water and liquid Nitrogen* Additionally, the physical properties of cryogenic fluids, other than vapor pressure, are also thermosensible, as illustrated in Figure 3. Thus, from a standpoint of numerical computations, simulating cryogenic cavitation implies a tight coupling between the non linear energy equation, momentum equations, and the cavitation model, via the iterative update of fluid properties (such as vapor pressure, densities, specific heat, thermal conductivity, viscosity etc.) with changes in the local temperature. Encountering these difficulties is expected to yield a numerical methodology specifically wellsuited for cryogenic cavitation, and forms the key emphasis of the present study. Obtained from REFPROP v 7.0 by Lemmon et al. (2002) 1 i Tomi ~ u i~ K) Trmp~rfrtur (K) Tmparalur (K) Twmprtaluf (K) 1.3 Contributions of the Current Study The major purpose of the present study is to develop a robust and comprehensive computational tool to simulate cavitating flow under cryogenic conditions. The specific contributions of the endeavor are summarized as follows: (a) A review of the experimental and computational studies on cryogenic cavitation (b) Coupling of energy equation to the existing cavitation framework in conjunction with iterative update of the real fluid properties with respect to the local temperature (c) Adaptation of an existing cavitation model (Senocak and Shyy 2004a, 2004b) to accommodate the physics of the mushy nature of cavitation observed in cryogenic fluids (Hord 1973, Sarosdy and Acosta 1961) (d) Demonstration of the impact of thermodynamic effects on cavitation over wide ranging temperatures, for two different cryogenic fluids. Assessment of the computational framework alongside available experimental and numerical data. (e) Global sensitivity analysis of the computational predictions (pressure and temperature) with respect to the cavitation model parameters and the temperature dependent material properties, via employing the response surface approach. (f) Adaptation of the pressurebased operator splitting method, PISO (Issa 1985), to multiphase environments typically characterized by strong interactions between the governing equations and steep variations of flow variables across the phase boundary (g) Assessment of the stability and accuracy of the noniterative algorithm (PISO variant) on the test problem of Gallium fusion. (h) Timedependent computations of cryogenic cavitation (with the PISO variant) by applying perturbation to the inlet temperature (i) Employment of Proper Orthogonal Decomposition (POD) to offer a concise representation to the simulated CFD data. CHAPTER 2 LITERATURE REVIEW Cavitation has been the focal point of numerous experimental and numerical studies in the area of fluid dynamics. A review of these studies is presented in this chapter. Since cryogenic cavitation remains the primary interest of this study, the general review on cavitation studies is purely restricted according to relevance. Specifically, the numerical approaches in terms of cavitation, compressibility, and turbulence modeling are briefly reviewed in the earlier section with reference to pertinent experiments. The later section mainly delves into the issues of thermal effects of cavitation. Current status of numerical strategies in modeling cryogenic cavitation, and their merits and limitations are reported to underscore the gap bridged by the current research study. 2.1 General Review of Recent Studies Computational modeling of cavitation has complemented experimental research on this topic for a long time. Some earlier studies (Reboud et al. 1990, Deshpande et al. 1994) relied on potential flow assumption (Euler equations) to simulate flow around the cavitating body. However, simulation strategies by solving the NavierStokes equations have gained momentum only in the last decade. Studies in this regard can be broadly classified based on their interface capturing method. Chen and Heister (1996) and Deshpande et al. (1997) adopted the interface tracking Marker and Cell approach in their respective studies, which were characterized by timewise grid regeneration and the constant cavitypressure assumption. The liquidvapor interface in these studies was explicitly updated at each time step by monitoring the surface pressure, followed by its reattachment with an appropriate wake model. This effort was mainly wellsuited for only sheet cavitation. The second category, which is the homogeneous flow model, has been a more popular approach, wherein the modeling for both phases is adopted via a singlefluid approach. The density change over the interface is simply modeled by a liquid mass fraction (f) or a liquid volume fraction (a,) that assumes values between 0 and 1. The mixture density can be expressed in terms of either fraction as follows: Pm =ap, +(1a)p, (2.1) = P + (2.2) The precise role of various cavitation models, which are reviewed later, is prediction of this volume/mass fraction as a function of space and time. Both, densitybased and pressurebased methods have been successfully adopted in conjunction with the singlefluid method in numerous studies. Due to unfeasibility of LES or DNS methods for multiphase flow, RANS approach through ke turbulence model has been mostly employed in the past studies. The main limitation of densitybased methods (Merkle et al. 1998, Ahuja et al. 2001, Lindau et al. 2002, Iga et al. 2003) is requirement of preconditioning (Kunz et al. 2000) or the artificial density approach for flows which may be largely incompressible. Pressurebased methods (Ventikos and Tzabiras 2000, Athavale et al. 2001, Senocak and Shyy 2002, Singhal et al. 2002) on the other hand are applicable over a wide range of Mach numbers. Modeling the speed of sound in the mixture region and computational efficiency for unsteady calculations are the main issues for pressurebased solvers. A broad classification of the numerical methods is illustrated in Figure 4. Interface Capturing Strategy Explicit Tracking Methods Homogenous Flow (Marker and Cell) Model Euler Equations Navier Stokes (Potential Flow Assumption) Solver Artificial Compressibility Based or Pressure Based Methods Figure 4. General classification of numerical methods in cavitation Within the broad framework of above methods, several cavitation, compressibility, and turbulence models have been deployed, which are reviewed next as per the aforesaid order. As mentioned earlier, the role of cavitation modeling is basically determination of the phase fraction. Different ideas have been proposed to generate the variable density field. Some studies solved the energy equation and obtained the density either by employing an equation of state (EOS) (Delannoy and Kueny 1990, Edwards et al. 2000) or from thermodynamic tables (Ventikos and Tzabiras 2000). Saurel et al. 1999 developed a cavitation model for hypervelocity underwater projectiles based on separate EOS for the liquid, vapor and mixture zones. EOS approach in first place does not capture the essential cavitation dynamics due to its equilibrium assumption in the phase change process. Furthermore, prevalence of isothermal condition in case of non thermosensible fluids such as water imparts a barotropic form, Pm = f(p), to the EOS. Thus, cavitation modeling via employing the EOS is devoid of the capability to capture barotropic vorticity generation in the wake region as demonstrated by Gopalan and Katz (2000). Transportequation based cavitation models in contrast overcome the above limitations and are more popular in research studies. Typically, this approach determines the liquid volume fraction (a,) or the vapor mass fraction (f,) by solving its transport equation as shown below. S+ V.(cai) = i+ + i (2.3) at a + V.(p, fii) = ih + (2.4) at Formulation of the source terms shown in above equation(s) constitutes the major effort in model development. Singhal et al. (1997), Merkle et al. (1998), Kunz et al. (2000), and Singhal et al. (2002) formulated these source terms strongly based on empirical judgment. However, Senocak and Shyy (2002, 2004a) developed a cavitation model fundamentally relying on interfacial mass and momentum transfer. Though their model was not completely empiricismfree, it transformed the empirical coefficients used in the earlier models into a physically explicable form. The ability of the model to capture the 1 barotropic vorticity, Vp xV in the closure region was also clearly demonstrated by P Senocak and Shyy (2002, 2004a). The source terms of each of the above models along with value of empirical constants are tabulated in Table 2. Table 2 Source terms in cavitation models Authors m' Production Term m Destruction term Singhal et al.(1997) CprodAX(p p,O0)(l a,) CdspMIN(p p,,0)a, (0.5p,U~)t. (0.5p,/.)pt Merkle et al. (1998) Cr = 8x101 Cds =1 Kunz et al. (2000) Cprodc a (1 a,) CppMIN(p p,,0)a, pit (0.5pU )Pit, Cprod = 3x 104 Cd, =1 Singhal et al.(2002) U 2 2AX( p,,0) 1/2 U 2MIN(p p, 0)1/2 pr 1 v dest r Iv 7 3 PI 3 Cprod = 3.675 x 103 Cds / = 1.225 x 103 Senocak and Shyy CprodAX (p p,0)(1 al) CdtpMIN(p p,,0)a, (0.5pU )t (0.5P/U )pt (2002, 2004a) Cro 1 Cd 1 0.5p,U, (P )U ,)2 0.5pU (P PV)(U,, U,)2 As clearly seen from the above table, the cavitation model of Senocak and Shyy (2002) is consistent with that of Singhal et al. (1997) and Merkle et al. (1998). However, the model constants have assumed physicality. The terms U,,, & U, ,, which represent the liquid vapor interface velocity and the normal component of interfacial vapor velocity are calculated by suitable approximations (Senocak and Shyy 2002, 2004a). Numerous experimental studies (Leighton et al. 1990, Clarke and Leighton 2000) addressing compressible bubble oscillations are available. Muzio et al. (1998) developed a numerical compressibility model for acoustic cavitation in a bubble inclusive of viscosity and surface tension effects. However, these models are difficult to deploy for practical multiphase computations. Figure 5 (adapted from Hosangadi et al. 2002) illustrates the modeled behavior of speed of sound v/s the phase fraction. As observed from the figure, the speed of sound in the biphasic mixture can be 23 orders of magnitude lower than the individual phases. Pure liquid 1400 1200 S1000 800 I 600  400 Pure vapor 200  0 0.25 0.5 0.75 1 Gas fraction Figure 5. Variation of Speed of Sound with phase fraction Consequently, a bulk incompressible flow may transform into a transonic or supersonic stage locally in the mixture region. Due to lack of dependable equation of state for multiphase mixtures, modeling sound propagation, which is an imperative issue in the numerical computations, is still an open question. A closed form expression for the speed of sound in the mixture region may be obtained by eigenvalue analysis on the strongly conservative form of the governing equations (Hosangadi et al. 2002), as shown below. m= P a + a2] (2.5) Cm pva pa, Venkateswaran et al. (2002) proposed use of perturbation theory to obtain an efficient preconditioned system of equations, which were consistent in the incompressible as well as compressible regime. Improvement in cavitation dynamics by accounting the compressibility effects was reported. The incorporation of Speed of Sound (SoS) model into pressurebased cavitation computations was significantly advanced by Senocak and Shyy (2002, 2003). In pressurebased solvers, the SoS model affects the solution mainly through the pressure correction equation. The following relationship is adopted between the density correction and the pressure correction terms, while enforcing the mass conservation treatment through the pressure correction equation: P'm = Cp' (2.6) The implementation of equation (2.6) imparts a convectivediffusive form to the pressure equation. Two SoS models were proposed by Senocak and Shyy (2003, 2004a, 2004b). SoS1: C = =C(1 ) (2.7) SoS2: C = ( ) p) p1 (2.8) Ap Op P, /Pr, While SoS1 is a suitable approximation to the curve shown in Figure 5, SoS2 approximates the fundamental definition of speed of sound by adopting a central difference spatial derivative along the streamline direction (0 instead of differentiating along the isentropic curve. Computations by Senocak and Shyy (2003) on convergent divergent nozzle demonstrated far better capability of SoS2 to mimic the transient behavior observed in experiments. Wu et al. (2003b) extended the model assessment by pointing out the dramatic time scale differences between the two models. This points the fact that compressibility modeling is a sensitivity issue and must be handled carefully. RANSbased approaches in form of twoequation k e models have been actively pursued to model turbulent cavitating flows. The k and E transport equations along with the definition of turbulent viscosity are summarized below. a(pk) (piuk) [ u, tk + =P, PE+ [(+ ) at x, xa 0 xa a(P (P ) + E2 a U au + = c,, P PtCp +[(p+ t ] at ax, k k [Ox x] The production of turbulent kinetic energy (Pt) is defined as: p = 1R a' Pr=r while, the turbulent viscosity is defined as: t = p "k 2 (2.9) (2.10) (2.11) (2.12) The model coefficients, namely, C,1,C, ok, and oa have three known nontrivial variants which have been summarized in Table 3. Table 3. Variants of the k ; model Authors Cg, and Relevant Details C2, 0k C, Launder and 1.44 1.92 1.3 1.0 Spalding (1974) Shyy et al. (1997) 115 + 0.25 P, 1.9 1.15 0.89 Younis (2003) PQ k )( 1.9 1.15 0.89 u(1.15+0.25 )x(1+0.38 I/Q) (personal where Q = k + (u2 +v2 + w2)/2 communication) Johansen et al. 1.44 1.92 1.3 1.0 (2004) p C,, k2 A F, C F,C = 0.09,F = Min[1, C  s k While the Launder and Spalding (1974) model is calibrated for equilibrium shear flows, the model by Shyy et al. (1997) accommodates nonequilibrium effects by introducing a subtle turbulent time scale into C,1. The RANS model of Younis (2003, personal communication), in comparison, accounts for the time history effects of the flow. Wu et al. (2003a, 2003b) assessed the above RANS models on turbulent cavitating flow in a valve. Experimental visuals of the valve flow (Wang 1999) have demonstrated cavitation instability in form of periodic cavity detachment and shedding. However, Wu et al. (2003b) reported that the k e model over predicts the turbulent viscosity and damps such instabilities. Consequently, the RANS computations were unable to capture the shedding phenomenon, and also showed restrained sensitivity to the above variants of the k E model. From standpoint of an alternate approach, the impact of filterbased turbulence modeling on cavitating flow around a hydrofoil was reported (Wu et al. 2003c). The filterbased model relies on the twoequation formulation and uses identical coefficient values as proposed by Launder and Spalding (1974), but imposes a filter on the turbulent viscosity as seen below. pC, k2 pt = F, C = 0.09 (2.13) The filter function (F1 is defined in terms of filter size (A) as: Ac F = Min[1,CA ]; CA = 1 (2.14) Note that the forms of viscosity in equations (2.12) and (2.13) are comparable barring the filter function. The proposed model recovers the Launder and Spalding (1974) model for coarse filter sizes. Furthermore, at nearwall regions, the imposed filter value F = 1 enables the use of wall functions to model the shear layer. However, in the far field zone if the filter size is able address the turbulent length scale , the solution is computed s directly (/u, = 0 ). The filterbased model is also characterized by the independence of the filter size from the grid size. This model enhanced the prediction of flow structure for singlephase flow across a solid cylinder (Johansen et al. 2004). Wu et al. (2003c) reported substantial unsteady characteristics for cavitating flow around a hydrofoil because of this newly developed model. Furthermore, the timeaveraged results of Wu et al. (2003c) (surface pressure, cavity morphology, lift, drag etc.) are consistent to those obtained by alternate studies (Kunz et al. 2003, CoutierDelgosha et al. 2003, Qin et al. 2003). Further examination in the context of filterbased modeling and cavitating flows was performed over the ClarkY aerofoil and a convergentdivergent nozzle (Wu et al. 2004, 2005). The filterbased model produced pronounced timedependent behavior in either case due to significantly low levels of eddy viscosity. While the timeaveraged results showed consistency to experimental data, they were unable to capture the essence of unsteady phenomena in the flowfield such as wave propagation. In addition to implementing the twoequation model for turbulent viscosity, Athavale et al. (2000) also accounted for turbulent pressure fluctuations. Thus, the threshold cavitation pressure (P,) was modified as: P' = P, + /2 (2.15) The turbulent pressure fluctuations were modeled as follows: p, = 0.39pk (2.16) Though their computations produced consistent results, the precise effect of incorporating the turbulent pressure fluctuations was not discerned. In addition to the above review, Wang et al. (2001), Senocak and Shyy (2002, 2004a, 2004b), Ahuja et al. (2001), Venkateswaran et al. (2002), and Preston et al. (2001) have also reviewed the recent efforts made in computational and modeling aspects. 2.2 Modeling Thermal Effects of Cavitation Majority of studies on cavitation have made the assumption of isothermal conditions since they focused on water. However, as explained earlier, these assumptions are not suitable under cryogenic conditions because of their low liquidvapor density ratios, low thermal conductivities, and steep slope of pressuretemperature saturation curves. Efforts on experimental and numerical investigation of cryogenic cavitation are dated as back as 1969. Though the number of experimental studies on this front is restricted due to the low temperature conditions, sufficient benchmark data for the purpose of numerical validation is available. However, there is a dearth of robust numerical techniques to tackle this problem numerically. The following subsections provide fundamental insight into the phenomenon in addition to a literature review. 2.2.1 Scaling Laws Similarity of cavitation dynamics is dictated primarily by the cavitation number (o) defined as (Brennen 1994, 1995): a P ZPT) (2.17) 0.5 p,U Under cryogenic conditions, however, cavitation occurs at the local vapor pressure dominated by the temperature depression. Thus, the cavitation number for cryogenic fluids is modified as: PP,(T) (2.18) S 5p(2.18) 0.5pU where, Tc is the local temperature in the cavity. The two cavitation numbers can be related by a firstorder approximation as follows: 1 pU (o )= dPv(T T.) (2.19) 2 dT Clearly, the local temperature depression (T T.) causes an increase in the effective cavitation number, consequently reducing the cavitation intensity. Furthermore, equation (2.19) underscores the effect of the steep pressuretemperature curves shown in Chapter 1. Quantification of the temperature drop in cryogenic cavitation has been traditionally assessed in terms of a nondimensional temperature drop termed as Bfactor (Ruggeri and Moore 1969). A simple heat balance between the two phases can estimate the scale of temperature difference caused by thermal effect. pt)L = pltlCPIAT (2.20) Here, vt and vu are volume flow rates for the vapor and liquid phase respectively. The B factor can then be estimated as: B= ;AAT pL (2.21) v, AT* P1CP Consider the following two flow scenarios for estimating B, as shown in Figure 6 (adapted from Franc et al. 2003). U0 (a) Twophase cavity U0 (b) 100% vapor cavity Thermal boundary Figure 6. Two cavitation cases for Bfactor analysis The estimation ofBfactor for the case (a) (twophase cavity) in above figure is expressed in following equation. a v0 a, U.9; v, (1 a,)U. 9; B = 1 (2.22) This points the fact that except for the pure vapor region B has an 0(1) value. For the case (b) in Figure 6 (adapted from Franc et al. 2003), where the cavity is assumed to be filled with 100% vapor, Fruman et al. (1991) provided an estimate of B based on thermal boundary layer effect as: B = (2.23) gT raL SU0 Here, a is the thermal or eddy diffusivity, Le is the cavity length, and '5 and ,5 represent the thickness of the cavity and the thermal boundary layer, respectively. It is evident from equation (2.23) that the temperature depression is also strongly dependent on thermal diffusivity and flow properties. The temperature scale (AT*) in case of water and LH2 has a value of 0.01 and 1.2 K, respectively (Franc et al. 2003). The difference in these values provides an assessment of the pronounced thermal effects in LH2. Thus, by knowing the values of AT* and B (from equations (2.22) and (2.23)), the actual temperature drop can be estimated. The Bfactor, however, fails to consider time dependent or transient thermal effects due to its dependence on a steady heat balance equation. Furthermore, the sensitivity of vapor pressure to the temperature drop, which is largely responsible in altering cavity morphology (Deshpande et al. 1997), is not accounted by it. As a result, though the Bfactor may estimate the temperature drop reasonably, it is inadequate to evaluate the impact of the thermal drop on the cavity structure and the overall flow. Brennen (1994, 1995) developed a more appropriate parameter to assess the thermodynamic effect by incorporating it into the RayleighPlesset equation (equation(2.24)) for bubble dynamics. d2R 3 dR P1[R + 2 (R) = pv(T) p. (2.24) dt 2 dt With help of equation(2.19), we can rewrite above equation as: d2R 3 dR dp p [R d + 3 )2]+ d AT=p(To) p (2.25) dt 2 dt dT From standpoint of a transiently evolving bubble, the heat flux q at any time t can be expressed as: AT q=K K (2.26) VaZ The denominator in equation(2.26), vt, represents the thickness of the evolving thermal boundary layer at time t. The heat balance across the bubble interface is expressed as follows. d4 q4rR2 = pL d[ 4 rR3] (2.27) dt 3 Combining equations (2.26) and (2.27) we obtain: AT Rv pL (2.28) 'Fa pCPl Introducing this temperature difference into equation(2.25) we obtain: d2R 3 dR2 P,(T.) P [Rd +2 ( ) 2]+ jXR (T= (2.29) dt2 2 dt pl A close observation yields the fact that the impact of thermal effect on bubble dynamics depends on Y which is defined as: Y= pL d= (2.30) p2Cp, a dT The units of are m/s3/2 and it proposes a criterion to determine if cavitation process is thermally controlled or not. Franc et al. (2003) have recently extended the above analysis to pose a criterion for dynamic similarity between two thermally controlled cavitating flows. They replaced the time dependency in equation(2.29) by spatialdependency through a simple transformation x = Ut Here, x is the distance traversed by a bubble in the flowfield in time t. If D is chosen as a characteristic length scale of the problem, the RayleighPlesset equation can be recast in the following form: 3 _V C +U [RR+3RR2]+ DR!= (2.31) 2 U 2 In equation (2.31), all the quantities with a bar are nondimensional, and Cp is the pressure coefficient. All the derivatives are with respect to x = x/D The above equation points out the fact that, besides o., two thermally dominated flows can be dynamically similar if they have a consistent value of the nondimensional quantity X D/UU It is important to note the important role of velocity scale U. in the quantification of thermal effect at this juncture. Franc et al. (2003) suggested that though thorough scaling laws for thermosensible cavitation are difficult to develop, a rough assessment may be gained from above equation. However, it is imperative to highlight that all the above scaling laws have been developed using either steadystate heat balance or single bubble dynamics. As a consequence, their applicability to general engineering environments and complex flow cases is questionable. The following sections will discuss the experimental and numerical investigations of cavitation with thermal consideration. Particularly, emphasis is laid on the limitations of currently known numerical techniques. Experimental studies are cited solely according to their relevance to the current study. 2.2.2 Experimental Studies Sarosdy and Acosta (1961) detected significant difference between water cavitation and Freon cavitation. Their apparatus comprised a hydraulic loop with an investigative window. While water cavitation was clear and more intense, they reported that cavitation in Freon, under similar conditions, was frothy with greater entrainment rates and lower intensity. Though their observations clearly unveiled the dominance of thermal effects in Freon, they were not corroborated with physical understanding or numerical data. The thermodynamic effects in cavitation were experimentally quantified as early as 1969. Ruggeri et al. (1969) investigated methods to predict performance of pumps under cavitating conditions for different temperatures, fluids, and operating conditions. Typically, strategies to predict the Net Positive Suction Head (NPSH) were developed. The slope of pressuretemperature saturation curve was approximated by the Clausius Clapeyron equation. Pump performance under various flow conditions such as discharge coefficient and impellor frequency was assessed for variety of fluids such as water, LH2, and butane. Hord (1973a, 1973b) published comprehensive experimental data on cryogenic cavitation in ogives and hydrofoils. These geometries were mounted inside a tunnel with a glass window to capture visuals of the cavitation zone. Pressure and temperature were measured at five probe location over the geometries. Several experiments were performed under varying inlet conditions and their results were documented along with the instrumentation error. As a result, Hord's data are considered benchmark results for validating numerical techniques for thermodynamic effects in cavitation. From standpoint of latest investigations, Fruman et al. (1991) proposed that thermal effect of cavitation can be estimated by attributing a rough wall behavior to the cavity interface. Thus, heat transfer equations for a boundary layer flow over a flat plate were applied to the problem. The volume flow rate in the cavity was estimated by producing an airventilated cavity of similar shape and size. An intrinsic limitation of the above method is its applicability to only sheettype cavitation. Larrarte et al. (1995) used highspeed photography and video imaging to observe natural as well as ventilated cavities on a hydrofoil. They also examined the effect of buoyancy on interfacial stability by conducting experiments at negative and positive AOA. They reported that vapor production rate for a growing cavity may differ substantially from the vapor production rate of a steady cavity. Furthermore, there may not be any vapor production during the detachment stage of the cavity. They also noticed that cavity interface under the effect of gravity may demonstrate greater stability. Fruman et al. (1999) investigated cavitation in R114 on a venturi section employing assumptions similar to their previous work (Fruman et al. 1991). The temperature on the cavity surface was estimated using the following flat plate equation. q 2.1 Te = T+ [1 2 (1 Pr)] (2.32) ple 0.5plCUC Re X Note that Cf is the coefficient of friction and Tpat is assumed to be equal to the local cavity surface temperature. The heat flux q on cavity surface was estimated as: q =pLUoC, (2.33) The discharge coefficient C. was obtained from an airventilated cavity of similar size. They estimated the temperature drop via the flat plate equation and measured it experimentally as well. These two corresponding results showed reasonable agreement. Franc et al. (2001) employed pressure spectra to investigate R114 cavitation on inducer blades. The impact of thermodynamic effect was examined at three reference fluid temperatures. They reported a delay in the onset of blade cavitation at higher reference fluid temperatures, which was attributed to suppression of cavitation by thermal effects. Franc et al. (2003) further investigated thermal effects in a cavitating inducer. By employing pressure spectra they observed shift in nature of cavitation from alternate blade cavitation to rotating cavitation with decrease in cavitation number. The earlier is characterized by a frequency 2fb (fb is the rotor frequency with 4 blades) while the later is characterized by a resonant frequency fb. Furthermore, R114 was employed as the test fluid with the view of extending its results to predicting cavitating in LH2. The scaling analysis provided in section 2.2.1 was developed by Franc et al. (2003) mainly to ensure dynamic similarity of cavitation in their experiments. 2.2.3 Numerical Modeling of Thermal Effects Numerical modeling has been implemented in cavitation studies broadly for two thermodynamic aspects. Firstly, attempts to model the compressible/pressure work in bubble oscillations have been made. Lertnuwat et al. (2001) modeled bubble oscillations by applying thermodynamic considerations to the RayleighPlesset equation, and compared the solutions to full DNS calculations. The bubble model showed good agreement with the DNS results. The modeled behavior, however, deviated from the DNS solutions under i\,theli imal and adiabatic assumptions. bubble liquid Figure 7. Schematic of bubble model for extracting speed of sound Rachid (2003) developed a theoretical model for accounting compressive effects of a liquidvapor mixture. The actual behavior of the mixture along with the dissipative effects associated with phase transformation was found to lie between two limiting reversible cases. One in which phase change occurs under equilibrium at a constant pressure, and the other in which the vapor expands and contracts reversibly in the mixture without undergoing phase change. Rapposelli and Agostino (2003) recently extracted the speed of sound for various fluids such as water, LOX, LH2 etc. employing a bubble model and rigorous thermodynamic relationships. The control volume (V) of the bubble is illustrated in Figure 7 (adapted from Rapposelli and Agostino 2003), and can be expressed as: V = (1 El)V + (E)V + (1 E,)V + (,)V, (2.34) The model assumed that thermodynamic equilibrium between the two phases is only achieved amidst fractions E, and E, of the total volume. Subsequently, the remaining fractions of the two phases were assumed to behave isentropically. If mi and m, are masses associated with the respective phases, the differential volume change dV can be expressed as: dV (1 a, d)p, d p, V p, dp' p, dp ( ) a, dp, 8 dp )S, (2.35) P, dp p, dp +dm(av O) m, m1 A close observation of above equation yields that the modeling extremities of E1, = 0 and E1,, =1 also correspond to the thermodynamic extremities mentioned in abovementioned analyses by Lertnuwat et al. (2001) and Rachid (2003) (italicized in the above description). Finally, substitution of various thermodynamic relations into equation (2.35) yields a thermally consistent speed of sound in the medium. Rapposelli and Agostino (2003) reported that their developed model was able to capture most features of bubble dynamics reasonably well. The second thermodynamic aspect of cavitation, which also forms the focal point of this study, is the effect of latent heat transfer. The number of numerical studies, at least in open literature, in this regard is highly restricted. Reboud et al. (1990) proposed a partial cavitation model for cryogenic cavitation. The model comprised three steps, which were closely adapted for sheet cavitation, in form of an iterative loop. (a) Potential flow equations were utilized to compute the liquid flow field. The pressure distribution on the hydrofoil surface was imposed as a boundary condition based on the experimental data. Actually, this fact led to the model being called 'partial'. The interface was tracked explicitly based on the local pressure. The wake was represented by imposing a reattachment law. (b) The vapor flow inside the cavity was solved by parabolized Navier Stokes equations. The change in cavity thickness yielded the increase in vapor volume and thus the heat flux at each section. (c) The temperature drop over the cavity was evaluated with the following equation: aT q = K, O (2.36) The value of turbulent diffusivity Kt in the computations was arbitrarily chosen to yield best agreement to the experimental results, and q was calculated from step (b). The iterative implementation of steps (a) (c) yielded the appropriate cavity shape in conjunction with the thermal effect. Similar 3step approach was adopted by Delannoy (1993) to numerically reproduce the test data with R113 on a convergentdivergent tunnel section. The main drawback of both the above methods is their predictive capability is severely limited. This is mainly because these studies do not solve the energy equation and depend greatly on simplistic assumptions for calculating heat transfer rates. Deshpande et al. (1997) developed an improved methodology for cryogenic cavitation. A preconditioned densitybased formulation was employed along with adequate modeling assumptions for the vapor flow inside the cavity and the boundary conditions for temperature. The interface was captured with explicit tracking strategies. The temperature equation was solved only in the liquid domain by applying Neumann boundary conditions on the cavity surface. The temperature gradient on the cavity surface was derived from a local heat balance similar to equation(2.36). The bulk velocity inside the vapor cavity was assumed equal to the free stream velocity. Tokumasu et al. (2002, 2003) effectively enhanced the model of Deshpande et al. (1997) by improving the modeling of vapor flow inside the cavity. Despite the improvements in the original approach (Deshpande et al. 1997), it is important to underscore the limitation that both the above studies did not solve the energy equation inside the cavity region. Hosangadi and Ahuja (2003, 2005), and Hosangadi et al. (2003) recently reported numerical studies on cavitation using LOX, LH2, and liquid nitrogen. Their numerical approach was primarily densitybased. Their pressure and temperature predictions over a hydrofoil geometry (Hord 1973a) showed inconsistent agreement with the experimental data, especially (Hord 1973a) at the cavity closure region. Furthermore, Hosangadi and Ahuja (2005), who employed the Merkle et al. (1998) model in their computations, suggested significantly lower values of the cavitation model parameters for the cryogenic cases as compared to their previous calibrations (Ahuja et al. 2001) for noncryogenic fluids. 31 The above review (summarized in Table 4) points out the limited effort and the wide scope for improving numerical modeling of thermal effects in cryogenic cavitation. A broadly applicable and more robust numerical methodology is expected to be a significant asset to prediction and critical investigation of cryogenic cavitation. Table 4. Summary of studies on thermal effects in cavitation Author and Year Method (Experimental/Numerical) Main Findings Sarosdy and Acosta Experimental Freon cavitates with a suppressed intensity 1961 Cavitation in hydraulic loop as compared to water Ruggeri et al. Experimental Provided assessment of pump NPSH under 1969 Cavitation in pumps cryogenic conditions Hord Experimental Published comprehensive test data on 1973 ogives and hydrofoil Fruman et al. Experimental Estimated temperature drop using flat plate 1991 Natural and ventilated cavities boundary layer equations Larrarte et al. Experimental Vapor production for growing cavities in 1995 High speed photography cryogenic fluids may not be estimated by ventilated cavities Fruman et al. Experimental Flat plate equations to model sheet 1999 Venturi section and R114 cavitation can give good results under certain restrictions Franc et al. Experimental Reported delay in onset of cavitation due 2001 Cavitation in pump inducers to thermal effect Franc et al. Experimental Provided scaling analysis to ensure 2003 Cavitation in pump inducers dynamic similarity Reboud et al. Numerical Semi empirical numerical model 1990 Potential flow equations Did not solve energy equation Suitable for sheet cavitation Delannoy Numerical Semiempirical numerical model 1993 Potential flow equations Did not solve energy equation Suitable for sheet cavitation Deshpande et al. Numerical Density based formulation 1997 Explicit interface tracking Simplistic approximation for cavity vapor flow Did not solve energy equation in the vapor phase Tokumasu et al. Similar to Deshpande et al. (1997) Improved flow model for cavity 2002, 2003 Applicable to only sheet cavitation Did not solve energy equation in the vapor phase Lertnuwat et al. Numerical Incorporated energy balance for bubble 2001 Model for bubble oscillations oscillations Good agreement with DNS simulation Rachid Theoretical Dissipative effects in phase 2003 Compression model liquid vapor transformation intermediate between two mixture extreme reversible thermodynamic phenomena Rapposelli & Agostino Numerical Employed thermodynamic relations to 2003 Model for bubble oscillations extract speed of sound for various liquids Hosangadi & Ahuja Numerical Solved energy equation in the entire 2003, 2005 Densitybased approach domain with dynamic update of material properties Some inconsistency with experimental results noted Significant change in the cavitation model parameters between noncryogenic and cryogenic conditions CHAPTER 3 STEADY STATE COMPUTATIONS This chapter firstly delineates the governing equations that are employed in obtaining steadystate solutions to various cases on cryogenic cavitation. Theoretical formulation/derivation and computational implementation of various models, namely cavitation, turbulence, compressibility, and thermal modeling aspects, are further highlighted. The boundary conditions and steadystate results yielded by the computational procedure are discussed in detail following the description of the basic framework. 3.1 Governing Equations The set of governing equations for cryogenic cavitation under the singlefluid modeling strategy comprises the conservative form of the Favreaveraged NavierStokes equations, the enthalpy equation, the k twoequation turbulence closure, and a transport equation for the liquid volume fraction. The masscontinuity, momentum, enthalpy, and cavitation model equations are given below: 8pm a(pi.,) a + =0 (3.1) at ax O(p.ur) 0P uzu) ap a au, au 2 al/ ( + [,) ((pU +)( + 3)] (3.2) at ax ax 9x 9x 9x( 3 (3 a[p (h+ fL)]+a[pm, (h+ fL)]= [(R + ) ) (3.3) at 9x ix Pr Pr /x c a,+ (a + h (3.4) at dx We neglect the effects of compressible work and viscous dissipation from the energy equation because the temperature field in cryogenic cavitation is mainly dictated by the phenomenon of evaporative cooling. Here, the mixture density, sensible enthalpy, and the vapor mass fraction are respectively expressed as: Pm = Pla + (la,) (3.5) h=CpT (3.6) f= (1 (3.7) Pm The general framework of the NavierStokes solver (Senocak and Shyy 2002, Thakur et al. 2002) employs a pressurebased algorithm and the finitevolume approach. The governing equations are solved on multiblock, structured, curvilinear grids. The viscous terms are discretized by secondorder accurate central differencing while the convective terms are approximated by the secondorder accurate Controlled Variations Scheme (CVS) (Shyy and Thakur 1994). The use of CVS scheme prevents oscillations under shocklike expansion caused by the evaporation source term in the cavitation model, while retaining second order of formal accuracy. Steadystate computations are performed by discounting the timederivative terms in the governing equations and relaxing each equation to ensure a stable convergence to a steady state. The pressurevelocity coupling is implemented through the SIMPLEC (Versteeg and Malalasekera 1995) type of algorithm, cast in a combined Cartesian contravariant formulation (Thakur et al. 2002) for the dependent and flux variables, respectively, followed by adequate relaxation for each governing equation, to obtain steadystate results. The physical properties are updated, as explained earlier, after every iteration. 3.1.1. Cavitation Modeling Physically, the cavitation process is governed by thermodynamics and kinetics of the phase change process. The liquidvapor conversion associated with the cavitation process is modeled through &' and i terms in Eq. (3.4), which respectively represent, condensation and evaporation. The particular form of these phase transformation rates, which in case of cryogenic fluids also dictates the heat transfer process, forms the basis of the cavitation model. Given below are the three modeling approaches probed in the present study. 3.1.1.1 Merkle et al. Model The liquidvapor condensation rates for this particular model are given as (Merkle et al. 1998): CdestplMin(O, p p,)a p,(O.5pUU )t SCprodMax(0, p p,)(1 a)) (0.5p,U )t. Here, Cdet =1.0 and Cprod =80.0 are empirical constants tuned by validating the numerical results with experimental data. The time scale (t.) in the equation is defined as the ratio of the characteristic length scale to the reference velocity scale (to = D/U ). 3.1.1.2 Sharp Interfacial Dynamics Model (IDM) The source terms of this model are derived by applying mass and momentum balance across the cavity interface and appropriately eliminating the unquantifiable terms (Senocak and Shyy 2002, 2004a, 2004b). u (a) (b) p I Figure 8. Schematic of cavity models (a) Distinct interface with vaporous cavity (Sharp IDM) (b) Smudged interface with mushy cavity (Mushy IDM) A schematic of the cavity model is illustrated in Figure 8(a). The model relies on the assumption of a distinctly vaporous cavity with a thin biphasic zone separating it from the pure liquid region. The physical form of the source terms is given below: lh pMin(0, p p, )a, P,(U,,n U1,n)2(pI P)()t9 Max(0, p p,)(l a,) (U, ,n U1,n)2 (pI Pv )t. Here, the normal component of the velocity is calculated as (Senocak and Shyy 2002, 2004a): Va' U =un; n= l (3.10) Vac, Due to implicit tracking of the interface, the interfacial velocity, U,,,, in unsteady computations needs modeling efforts. Previous studies simplistically expressed the interfacial velocity (U1,,) in terms of the vapor normal velocity (Senocak and Shyy 2002, 2004a) (U,,n). Alternate methods of modeling U1,, are discussed in the Chapter 4 in the context of timedependent simulations. In comparison, for the steady computations in this chapter, we impose U, = 0 . 3.1.1.3 Mushy Interfacial Dynamics Model (IDM) Experimental visualizations of cryogenic cavitation (Sarosdy and Acosta 1961; Hord 1973a, 1973b) have clearly indicated a mushy nature of the cavity. This salient characteristic of cryogenic cavitation solicits an adaptation of the existing cavitation model to reflect the same. We choose the Sharp IDM (Senocak and Shyy 2004a) in our analysis because of its stronger physical reasoning. The following discussion serves to revise the above cavitation model by reexamining its derivation and assumptions, to appropriately accommodate the features of cryogenic cavitation. Although we refer the reader to literature (Senocak and Shyy 2004a) for the detailed derivation of the earlier model, we underscore the differences between the two approaches in the proceeding description. Figure 8(b) representing Mushy IDM depicts a cavity where the vapor pressure is a function of the local temperature, and the nature of cavitation demonstrates a weak intensity and consequently less probability for existence of pure vapor phase inside the cavitation zone. In comparison, Figure 8(a) representing Sharp IDM depicts a cavity typically produced under regular conditions (no temperature effects), characterized by a thin biphasic region separating the two phases. We initiate our approach similar to Senocak and Shyy (2004a) by formulating the mass and momentum balance condition at the cavity interface, which is assumed to separate the liquid and mixture regions. We neglect the viscous terms and surface tension effects from equation (3.12) assuming high Re flows and large cavity sizes respectively. p (U,, (U,) = P.(U.," VU,) (3.11) P P, = P,(Un, V J,) p,(Up,. V,)2 (3.12) It is noteworthy that Senocak and Shyy (2004a) formulated the above equations midst the pure vapor and mixture regions while hypothesizing an interface between the two. The two models differ in terms of their interpretation of the cavity attributes, namely, the density and the velocity field within the cavitation zone, while adhering to the fundamental idea of mass/momentum balance. Since we attribute a "frothy" nature to the cavity, the term fi.n is largely expected to represent the mixture normal velocity. Accounting for this fact, we eliminate the U,,, term from equation (3.12) using equation (3.11) as shown below. P1 P = (Um, U,,) Pm (Um,, U,,)2 (3.13) P1 Further rearrangement progressively yields the following equations: P1 = ,(Um,, U,)2[1 ] p(U,, U,,)2 1 a, (3.14) a, P(P, P') P(P P') a = =(a, + a,) (3.15) P, (P, p)(U,,, U,1)2 p, (P p,)(U U1)2 p,(p, P)a, p+ (p, p)(1 a1) a = + (3.16) p, (P P)(U,,, U~,,)2 p,(P Pv)(Um ,, U1),n) From standpoint of formulating source terms for the a, transport equation, we firstly normalize the above equation by the overall convective timescale (t = D/U,). Secondly, we apply conditional statements on the pressure terms to invoke either evaporation or condensation depending on the local pressure and vapor pressure (refer to Senocak and Shyy 2002 & 2004a). Lastly, we assume that the volume rates for the individual phases are interchangeable barring the sign convention. For instance, the evaporation term, ih in the a, equation or a, equation would bear the same magnitude, but negative or positive sign, respectively. We thus obtain the following source terms from the above analysis: pMin(0, p p,)a Pm (UIm, U],n)Z(P, P )t (3.17) p Max(O, p p,)(1 a,) p, (U,,, U_~,)2 (pi p,)t The normal component of the mixture velocity is expressed as shown in equation (3.10) and consistent with the normal component of the vapor velocity employed by the Sharp IDM. A comparison with equation (3.9) shows that a factor of p/Pm, weakens the evaporation term while a factor of p, / p strengthens the condensation term of the Sharp IDM to yield the Mushy IDM. For instance, at a nominal density ratio of 100 and a, = 0.5, the value of p / pv and p / p is 1.82 10 2 and 1.82, respectively. In terms of order of magnitude, the difference between the evaporation terms of the two models is relatively more pronounced than the condensation terms, especially for high liquidvapor density ratios. Given this fact, we attempt to make the Mushy IDM consistent with the Sharp IDM by providing an exponential transition of the evaporation source term from one model to the other as a function of the phase fraction (a, ). Specifically, we employ the evaporation term from equation (3.17) as we get closer to the liquid region, and that from equation (3.9) as we get farther from the liquid region. In comparison, the condensation term from equation (3.17), which is mainly expected to be active in the region outlining the 'vaporous' portion of the cavity, is unaltered. In summary, the Mushy IDM is eventually formulated as shown below. p Min(0, p p,)a, p Max(O, p p,)(1 a,) P (Um U (U,n)2(p p, )t p( (3.18) AL _L + (1. OPl)e(a,)p p pL p P P, P Pm Note that / is a free parameter, which regulates the switch between the evaporation terms and warrants calibration for different fluids. Its typical value could be 0(0.1). 4 5 "o 4 Figure 9. Behavior of p/pP and p,/P vs. at for the two models; pl/P = 100 and / = 0.09 Figure 9 contrasts the densityratio terms ( pl / p and pl /p) of the Sharp IDM and Mushy IDM for the chosen parameter values. 3.1.2 Turbulence Modeling For the system closure, the original twoequation turbulence model with wall functions is presented as follows (Thakur et al. 2002): G(pOk) *(p uk) \ , k + = Pt PmE + [(U + ) ] (3.19) dt &x x Ok Fi 5e90eairo ^Ipad/;/+v.a o hetomdl;p/p=0 n :?=0.0 O.~k O. u, k)a O(P(E) dPuE) E E2 at UOE + = P, C2Pm + [(u + ) ] (3.20) at cx k k cxj oxj The turbulence production (P,) and the Reynolds stress tensor is defined as: au, Pt = 7'j 0 =pMu'u' x (3.21) 2Pk au au Pu 'u' m lut ( a' + 3 x ax, The parameters for this model, namely, C, = 1.44, C,2 =1.92, o =1.3, o, =1.0 are adopted from the equilibrium shear flow calibration (Launder and Spalding 1974). The turbulent viscosity is defined as: pmC k2 pt ,,C = 0.09 (3.22) It should be noted that the turbulence closure and the eddy viscosity levels can affect the outcome of the simulated cavitation dynamics especially in case of unsteady simulations (as reviewed in Chapter 2). In this aspect, parallel efforts are being made in the context of filterbased turbulence modeling (Johansen et al. 2004; Wu et al. 2003c, 2004, 2005), which has shown to significantly increase the timedependency in cavitating flows. We do not explore these techniques in the interest of steadystate simulations, which disregard timedependent phenomena. 3.1.3 Speed of Sound (SoS) Modeling Due to lack of dependable equation of state for liquidvapor multiphase mixture, numerical modeling of sound propagation is still a topic of research. We refer the reader to past studies (Senocak and Shyy 2003, 2004a, 2004b, Wu et al. 2003b) for modeling options, their impact and issues, and just outline the currently employed SoS model below. SoS= C, = C(1a,) (3.23) The density correction term in the continuity equation is thus coupled to the pressure correction term as shown below. P = Cp (3.24) Senocak and Shyy (2002 & 2004a) suggested an 0(1) value for the constant C to expedite the convergence of the iterative computational algorithm. However, their recommendation is valid under normalized values for inlet velocity and liquid density. Since we employ dimensional form of equations for cryogenic fluids, we suggest an 0(1 / U ) value for C, which is consistent with the above suggestion in terms of the Mach number regime. The speed of sound affects the numerical calculation via the pressure correction equation by conditionally endowing it with a convectivediffusive form in the mixture region. In the pure liquid region, we recover the diffusive form of the pressure equation. 3.1.4 Thermal Modeling The thermal effects are mainly regulated by the evaporative cooling process, which is further manifested by the temperature dependence of physical properties and vapor pressure. 3.1.4.1 Fluid property update In the present study, we subject all the physical properties, namely, p1, Pv, ,v,, Cp, K, and L to temperature dependence. 43 S.......... Dsnallty kl m) Enthlpy kJkg) Figure 10. Pressuredensity and pressureenthalpy diagrams for liquid nitrogen in the liquidvapor saturation regime (Lemmon et al. 2002). Lines denote isotherms in Kelvin. As indicated by Figure 10, the physical properties are much stronger functions of temperature than pressure, and can fairly assume the respective values on the liquid vapor saturation curve at a given temperature. We update these properties from a NIST database (Lemmon et al. 2002) at the end of a computational iteration. Thus, we generate a lookup table of physical properties for a particular temperature range as a pre processing step. Subsequently, for any temperaturebased update, the table is searched by an efficient bisection algorithm (Press et al. 1992) and the required physical property is obtained by interpolating between the appropriate tabular entries. 3.1.4.2 Evaporative cooling effects The energy equation, (Eq.(3.3)), is recast into the following temperaturebased form, by separating the latent heat terms onto the righthandside. a[PmCpT]+ [pmu CpT]= P[CP(P + P ' et 9xg xo, Pr, Pr txe a a (3.25) [EP/(f;L)]+ [pu(f~L)]} v t soue/s energy source/sink term As seen from equation (3.25), the 'lumped' latent heat terms manifest as a nonlinear source term into the energy equation and physically represent the latent heat transfer rate. The spatial variation of thermodynamic properties and the evaporative cooling effect are intrinsically embedded into this transportbased source term. We calculate the source term by discretizing the associated derivatives in concert with the numerical schemes applied to the terms on the lefthandside of the equation. 3.1.5 Boundary Conditions The boundary conditions are implemented by stipulating the values of the velocity components (obtained from the experimental data), phase fraction, temperature, and turbulence quantities at the inlet. Furthermore, at the walls, pressure, phase fraction, and turbulence quantities are extrapolated, along with applying the noslip and adiabatic condition on the velocity and temperature, respectively. Pressure and other variables are extrapolated at the outlet boundaries, while enforcing global mass conservation by rectifying of the outlet velocity components. In addition, we also hold the pressure at the reference pressure point constantly at the reference value (specified by the experiments). This is simply done by adjusting the linear coefficients of the pressure correction equation at that point, to yield zero correction, at every iteration. Symbolically, this is achieved by substituting Ap = 1, Ap = BP = 0 into the following linear equation at the point of interest: Ap pp = ZAn + Bp (3.26) We observe that this adjustment imparts robustness and stability to the computation. 3.2 Results and Discussion In this section, we firstly observe the nuances of the Mushy IDM on noncryogenic cases. The purpose is to distill the impact solely caused due to the source term change between the Sharp and Mushy IDM, and contrast it against commonly referred experimental and computational results. Later, we extend our computations to the variable temperature environment of cryogenic fluids. Specifically, we observe the effect of the temperature field on the nature of cavitation under similar conditions. Furthermore, we assess the Mushy IDM with available experimental data (pressure and temperature) and compare it with the alternative cavitation models. In the process of calibrating the cavitation models for cryogenic fluids, we perform a global sensitivity analysis to evaluate the sensitivity of the prediction to changes in material properties and model parameters. We offer discussion about optimizing the model performance based on our sensitivity study. 3.2.1 Cavitation in Noncryogenic Fluids As mentioned above, we register first the impact of the Mushy IDM on non cryogenic cases. In this initial exercise, we consider the influence of the mushy formulation only to the liquid boundary layer encompassing the cavity. Symbolically, we implement the following formulation for the noncryogenic cases. lh pMin(0, p p,)al p (Um,n U,,,)2(P, Pv)t p Max(O, p p,)(1 a,) P.(Un, U,,)2 (P P)t (3.27) (3.27) ifa,< 0.99 P = else P PA P Pv P P, if a,<0.99 p 1 elseP P P+ P+ P, To prevent sudden discontinuity in the volume transfer rates, we perform a geometric smoothing operation over the source terms. Two flow configurations are discussed here, namely, cavitating flow over a hemispherical projectile (timeaveraged experimental data by Rouse and McNown (1948) at Re=1.36x105) and cavitating flow over the NACA66MOD hydrofoil (timeaveraged experimental data by Shen and Dimotakis (1989) at Re = 2x106). NOSLIP OUTLET hemispherical projectile INLET ( NOSLIP NOSLIP INLET OUTLET NOSLIP Figure 11. Illustration of the computational domains for hemispherical projectile and NACA66MOD hydrofoil (noncryogenic cases) hydrofoil (NOSLIP) The computational domains for the two geometries are depicted in Figure 11. We consider two grids with 158x66 and 292x121 points for the hemispherical geometry, which is axisymmetric. In case of the NACA66MOD hydrofoil, we defer to the judgment of a previous study (Senocak and Shyy 2004a), and employ the finer grid from that study. (a)  Mushy 18x66 gd (b) Sharp IDM ...... ushy IDM 292x121 grid  Merkle et al. Model 1 Mushy IDM 292 x121 grid Mushy DM 0 Exp. data 1 Exp. dushy ID U Exp. data 05 05 0 0 0 0 0 \ 05 05  0 1 2 3 4 0 1 2 3 4 s/D s/D Figure 12. Pressure coefficients over the hemispherical body (o = 0.4); D is the diameter of the hemispherical projectile. (a) Impact of grid refinement for Mushy IDM (b) Comparison between pressure coefficients of different models on the coarse grid We investigate the sensitivity of Mushy IDM to grid refinement for the two hemispherical body grids through the surface pressure plots in Figure 12. Note that the refinement factor in the vicinity of cavitating region is roughly 2.53. Figure 12(a) indicates a noticeable, though modest, effect of the grid refinement on the surface pressure. This can be mainly attributed to the griddependent geometric smoothing operation, which smears the mushy formulation over a larger portion in the coarser grid. Nonetheless, both the grids produce solutions that match the experimental data reasonably well. Furthermore, it is important to mention that the nearwall nodes of the coarser grid lie in the loglayer of the turbulent boundary layer enabling appropriate use of the wall function (Versteeg and Malalasekera 1995), while the fine grid is overrefined for the wall spacings to assume values in the suitable range. As pointed out by Senocak and Shyy (2004a), in view of the wall function treatment (Thakur et al. 2002; Versteeg and Malalasekera 1995), the presence of vapor in the cavity may substantially reduce the appropriate range of the wall node spacings. As a consequence, we employ the coarser grid for further computations. Figure 12(b) contrasts the performance of the Mushy IDM with the Merkle et al. Model (1998) and Sharp IDM (2004a) at c = 0.4. The pressure coefficient predicted by the three models illustrates noticeable variations in the condensation region of the cavity. Specifically, the Mushy IDM, as seen from Figure 12 (b), tends to produce a sharper recovery of the surface pressure in the cavity closure zone. Additionally, the Mushy IDM produces a small pressure dip (below the vapor pressure value) at the cavity outset due to the lower evaporation rate and higher condensation rate. (a) Liquid Volume Fraction 1 0 08 05 (b) Liquid Volume Fraction S04 028 07 05 (c) Liquid Volume Fraction 1 0 08 04 *02 Figure 13. Cavity shapes and flow structure for different cavitation models on hemispherical projectile (o = 0.4). (a) Merkle et al. Model (b) Sharp IDM (c) Mushy IDM 49 These features are reflected by the cavity sizes depicted in Figure 13. Note that the Mushy IDM not only shrinks the cavity length, in comparison with the Sharp IDM, but also impacts the flow structure in the recirculation zone of the cavity closure region. The above findings support our approach in providing an appropriate transition between the two models and a careful calibration to the parameter /7. (a) (b) (a) Sharp IDM  Sharp IDM  Merkle et al. Model  Merkle et al. Model MushylDM Mushy lDM 1 2 Exp. data 1 2 U Exp. data 08 k08  06 06  02 02 0 0 0 2 02  0 4 04 0 01 02 03 04 05 06 07 08 09 1 0 01 02 03 04 05 06 07 08 09 1 x/D x/D Figure 14. Pressure coefficient over the NACA66MOD hydrofoil at two different cavitation numbers; D is hydrofoil chord length. (a) a = 0.91 (b) a = 0.84 We further simulate cavitating flow around the NACA66MOD hydrofoil at two cavitation numbers, namely, 0.84 and 0.91, while maintaining the mushy formulation only for the cavity boundary. Again, consistent behavior in terms of surface pressure is induced by the Mushy IDM for this geometry, at both the cavitation numbers (Figure 14). In summary, the appropriate source term modulation imposed by the Mushy IDM tends to influence the prediction of surface pressure and cavity size in manners consistent with the experimental observation. The above assessment motivates the further implementation and enhancement of the Mushy IDM (with respect to /) on cryogenic flow cases. 3.2.2 Cavitation in Cryogenic Fluids In this section, we investigate the Mushy IDM for the cryogenic situation. We perform computations on two geometries experimentally investigated by Hord (1973a, 1973b), namely, a 2D quarter caliber hydrofoil (Hord 1973a) and an axisymmetric 0.357 inch ogive (Hord 1973b). These geometries were mounted inside suitably designed tunnels in the experimental setup. NOSLIP INLET OUTLET SYMMETRY hydrofoil surface (NOSLIP) NOSLIP INLET I I OUTLET SYMMETRY ogive surface (NOSLIP) Figure 15. Illustration of the computational domain accounting the tunnel for the hydrofoil (Hord 1973a) and 0.357inch ogive geometry (Hord 1973b) (cryogenic cases). Figure 15 illustrates the computational domains employed for these cases. Note that the figure shows only planar slices of the domains, which also model the respective tunnel shapes to account for the significant blockage effects. The mesh for the hydrofoil and the ogive geometry respectively comprises 320x70 and 340x70 points. The mesh distribution is chosen to facilitate adequate resolution of the cavitation zone. Furthermore, the nearwall resolution over all the noslip planes (cavitating geometries and tunnel walls) accounts for deployment of wall functions (Thakur et al. 2002; Versteeg and Malalasekera 1995). Table 5. Flow cases chosen for the hydrofoil geometry. Fluid: Case name: Inlet temperature: Freestream Re: Cav. No. (o.): Liq. N2 283B 77.65 K 4.7x106 1.73 Liq. N2 290C 83.06 K 9.1x106 1.70 Liq. N2 296B 88.54 K 1.1x107 1.61 Liq. H2 248C 20.46 K 1.8x107 1.60 Liq. H2 249D 20.70 K 2.0x107 1.57 Liq. H2 255C 22.20 K 2.5x107 1.49 Source: Hord (1973a) Table 6. Flow cases chosen for the ogive geometry. Fluid: Case name: Inlet temperature: Freestream Re: Cav. No. ( ,): Liq. N2 312D 83.00 K 9.0x106 0.46 Liq. N2 322E 88.56 K 1.2x 107 0.44 Liq. H2 349B 21.33 K 2.3x107 0.38 Source: Hord (1973b) Statisticallyaveraged pressure and temperature data are available for the two geometries at five probe locations over the body surfaces. The experimental findings report varying amounts of unsteady behavior in the cavity closure regions, although no casespecific information or data/visuals are available in that context. Hord conducted a series of experiments over both the geometries, using liquid nitrogen and hydrogen, by varying the inlet velocity, temperature, and pressure. In our computations, we have selected several cases, referenced alphanumerically in the reports by Hord (1973a, 1973b), with different freestream temperatures and cavitation numbers (see Table 5 and Table 6). (a) Computed pressure (b) Computed pressure 1 U Hord, 1973 (incipient data) 1 5 Hord, 1973 (incipient data) 05 1  0  05  05 05 25 1 5 I I I I I I 2 I I I 1 0 1 2 3 4 5 0 2 4 x/D x/D Figure 16. Noncavitating pressure distribution (a) case '290C', D represents hydrofoil thickness and x represents distance from the circular bend (b) case '312D', D represents ogive diameter and x represents distance from the leading edge To validate the use of real fluid properties, we obtain the singlephase flow solution to cases '290C' (hydrofoil; Re=9.1x106, c=1.7) and '312D' (ogive; Re=9.0x106, o,=0.46), and compare the computed surface pressure with the experimentally measured pressure under incipient (virtually noncavitating) conditions. Figure 16 demonstrates good agreement between the numerical and experimental data for both the geometries, and corroborates the correct input of physical properties. 3.2.2.1 Sensitivity analyses In general, cryogenic computations are prone to uncertainty due to a multitude of inputs, in contrast to the noncryogenic conditions. The cavitation model parameters, namely, Cdet, Cprod, and t, have been selected largely based on the noncryogenic conditions. Furthermore, the solutions can exhibit substantial sensitivity with respect to minor changes in the flow environment. For example, the uncertainties involved in the 53 temperaturedependent material properties may also cause noticeable differences in predictions. To address relevant issues in this context and confirm/improve our calibration of the Merkle et al. model constants, we perform a comprehensive Global Sensitivity Analysis (GSA) over the chosen case '290C'. We initiate our computations on the case '290C' (Re = 9.1 x106; o = 1.7) which is centrally located in the temperature range. We note that the previously calibrated values of the Merkle et al. Model (Cdet =1.0 and C prod=80.0) are inadequate to provide a good match with the experimental data under the cryogenic condition. This fact was also lately noted by Hosangadi and Ahuja (2005), who suggested lower values of cavitation model parameters in case of cryogenic fluids. C=0 068, C d=54 4544  C = 0, Cp 80 0 845 C = 68, Cprod Hord, 1973,p ... .. dest=l 0, Cprod80 0 30 84rd 1973 Hord, 1973, T 20 835  .  835 20 i 825. 82 0. 0' 81 5  81 10 805 0 5 0 05 1 15 0 05 1 15 2 x (inches from bend) x (inches from bend) Figure 17. Sensitivity of Merkle et al. Model prediction (surface pressure and temperature) to input parameters namely Cdet and Cprod for the hydrofoil geometry Consistently, in the present study, we elicit Cdet = 0.68 and Cprod = 54.4 via numerical experimentation, as more appropriate model parameters. Of course, such choices are empirically supported and need to be evaluated more systematically. To initiate such evaluations, Figure 17 portrays the response of the surface pressure and temperature to the revision of the model parameters. We choose Cde,, to, Pv, and L as design variables for the GSA, while holding the Re and ro constant for the given case. The chosen cavitation model parameters, namely Cdet and to, are perturbed on either side of their reference values (Cdt = 0.68; Cprod = 54.4) by 15%. In comparison, the material properties are perturbed within 10% of the value they assume from the NIST database (Lemmon et al. 2002), at every iteration. Given these ranges on the variables, we generate a design of experiments with 50 cases using a combination of Orthogonal Arrays (Owen 1992) and Face Centered Cubic Design (JMP 2002). RMS values of the hydrofoil surface pressure coefficient (Cp = (p pm)/(0.5p1U2)) and temperature, which are postprocessed from the CFD data, are selected as the objective functions. Subsequently, the two objectives are modeled by a reducedquadratic response surface through a leastsquares regression approach (JMP 2002). The coefficient of multiple regression (Myers and Montgomery 1995) in case of the pressure and temperature fit is 0.992 and 0.993 respectively, while the standard error is less than 1%. The fidelity of the response surfaces is also confirmed against 4 test data points. S 1.675 + 0.077C, 0.061p 0.082t: 0.007C*,dp  (3.28) 0.012C tst +0.011p*2 +0.009p9t* 0.004L2 +0.013t.2 s = 82.537 0.243Cd, + 0.107p* 0.048* + 0.220t + 0.024C*t 0.016Cd ,pv 0.015C,tL* +0.013pl L + (3.29) 0.018p tc +0.019L*t Here, CPM and Ts are the RMS values of the surface pressure coefficient and the surface temperature, respectively. The superscript, *, represents the normalized values of the design variables. Equations (3.28) and (3.29) represent the respective response surfaces expressed in terms of the normalized design variables. The surface coefficients upfront indicate the importance of the cavitation model parameters (Cde, and t.) and vapor density, and the insubstantial contribution of latent heat (L), to the variability in our objectives. We quantify these overall contributions by employing the variancebased, nonparametric global sensitivity method proposed by Sobol (1993). This method essentially comprises decomposition of the response surface into additive functions of increasing dimensionality. This allows the total variance in the data to be expressed as a combination of the main effect of each variable and its interactions with other variables (refer to Appendix A for a brief mathematical review). (a) (b) 42% m Latent heat Latent heat 4200 1 10 m Vapor density 46% m Vapor density oC dest oC dest Ot infinity 43% [ t infinity u% I' 38% Figure 18. Main contribution of each design variable to the sensitivity of Merkle et al. (1998) model prediction; case '290C' (a) Surface pressure (b) Surface temperature We implement the procedure expounded by Sobol (1993) on the response surfaces in equations (3.28) and (3.29) to yield the plots in Figure 18. The piecharts in the figure illustrate the percentage contribution of the main effect of each variable; since we find negligible variability due to the variable interactions. The charts firstly underscore the sensitivity of pressure and temperature predictions to the cavitation model parameters (Cde,, and t.). Secondly, the impact of p, is noticeable, while that of L is 56 insubstantial. These observations indicate that the design variables, unlike L, which appear either in mi or i+ may tend to register greater influence on the computed results. Thus, intuitively, U, and p,, which are omitted from the present GSA, are expected to induce large variability in the computation, as compared to other omitted properties such as K and Cp.  Merkle et al. Model * Hord,1973;p 15 0 05 1 x (inches from bend) Figure 19. Pressure and temperature prediction for Merkle et al. Model for the case with best match with experimental pressure; Cde, = 0.85; t = 0.85; p =1.1; L = 0.9  Merkle et al. Model Hord,1973;p 1  Merkle et al. Model Hord,1973;T I I II80I I I U U5 1 x (inches from bend) 15 2 Figure 20. Pressure and temperature prediction for Merkle et al. Model for the case with best match with experimental temperature; Cet =1.15; t = 0.85; p* =1.1;L* = 1.1 30 20 E 10 0 0. 10 20  Merkle et al. Model Hord,1973;T 0U U 05 1 x (inches from bend) 15 2 20.9.. Ub U Ub 1 x (inches from bend) I I I I I I I III II I I I Furthermore, the impact is expected to be consistent on pressure and temperature, as depicted by Figure 18, due to the tight coupling between various flow variables. However, it is important to mention that our observation is meant to illustrate the relative impact of several parameters. We also utilize the available data from our 50 cases to seek possible improvement of the Merkle at al. model parameters. Figure 19 illustrates the case which produces the least RMS error between the computed surface pressure and the experimental data. Conversely, Figure 20 portrays the case which produces the least RMS error between the computed surface temperature and the experimental data. These figures demonstrate that a single set of parameters may not provide optimal results for both pressure and temperature within the framework of current cavitation models and cryogenic conditions. This effort solicits multiobjective optimization strategies and deserves a separate study. However, we do report from the calculated error norms of the 50 cases that Cdet = 0.68 and Cprod = 54.4 provide the best balance between the temperature and pressure predictions for the chosen case. As a result, we hereafter employ these values for the Merkle et al. Model. We perform a simpler sensitivity analysis over the Mushy IDM since it has only one control parameter (/8). 58  = 0.07  =0.07  = 0.09  p=0.09  P= 0.11 ....... p=0.11 Hord,1973;p 84 5 Hord,1973;T 84  20  835  E./ / / 83 S82/ / 82/ / II82  81 5 81 80 5 20 80 05 0 05 1 15 0 05 1 15 x (inches from bend) x (inches from bend) Figure 21. Sensitivity of Mushy IDM prediction for case '290C' (surface pressure and temperature) to the exponential transitioning parameter 0 Figure 21 depicts the results obtained over various values of f/. We calibrate / = 0.09 from the observed results similarly based on a reasonable balance between pressure and temperature. Furthermore, Figure 21 suggests that the limiting case of = 0, which essentially recovers the evaporation term of the Sharp IDM, would substantially over predict the cavity size. This fact endorses the need to regulate the mass transfer rates modeled by the Sharp IDM, and subsequently the purported employment of the Mushy IDM. It is worthwhile to emphasize that the Mushy IDM manifests the regulation of mass transfer rates in cryogenic conditions incorporated empirically into the Merkle et al. Model (section 3.2.1; Ahuja et al. 2001; Hosangadi and Ahuja 2005) largely under physical pretext. 3.2.2.2 Assessment of cryogenic cavitation models over a wide range of conditions We further perform computations for all the other cases by employing both the Merkle et al.'s and the present Mushy IDM cavitation models.  Mushy DM  Merkle et al. Model U Hord, 1973;p 283B 5 0 05 1 x (inches from bend) 15  Mushy DM  Merkle et al. Model Hord, 1973;p A Hosangadi &Ahuja, 2005 290C 5 0 05 1 15 x (inches from bend)  MushylDM  Merkle et al. Model U Hord, 1973;p 296B 20 .... . 05 U 05 1 x (inches from bend) 1  Mushy DM  Merkle et al. Model U Hord,1973;T 283B 0 05 1 15 2 25 3 x (inches from bend)  Mushy DM  Merkle et al. Model U Hord, 1973; T A Hosangadi &Ahuja, 9 S290C A 290C 2005 ] I . I . I . I I I . I . I 0 05 1 15 2 25 3 x (inches from bend)  Mushy IDM  Merkle et al. Model Hord,1973;T 296B 0 05 1 15 2 25 3 x (inches from bend) Figure 22. Surface pressure and temperature for 2D hydrofoil for all cases involving liquid Nitrogen. The results referenced as 'Mushy IDM' and 'Merkle et al. Model' are contributions of the present study. I I IIII.. .I I Figure 22 contrasts the resultant surface pressure and temperature obtained with the hydrofoil geometry for the cases involving liquid Nitrogen. We firstly note that both the models are able to provide a reasonable balance between pressure and temperature from standpoint of their predictive capabilities. The differences between the experimental data and the two models are more pronounced than the mutual differences between the two models. Furthermore, the agreement with experimental data is better in case of pressure than in case of temperature, which is generally underpredicted at the leading probe point by both the models. It is also observed that both models produce a slight temperature rise above the reference fluid temperature at the cavity rear end, which is attributed to the release of latent heat during the condensation process. The Merkle et al. Model also produces a steeper recovery of pressure, as compared to the Mushy IDM, in the condensation region of the cavity, for the cases shown. This suggests higher/faster condensation rates for the Merkle et al. Model than the Mushy IDM. Secondly, we assess our results for the case '290C' along with latest computational data (Hosangadi and Ahuja 2005). Note that Hosangadi and Ahuja (2005) employed the Merkle et al. Model, adapted in terms of the vapor volume fraction (ca,), with substantially higher values of the model coefficients (Cdet =Cprod =100). The impact of these higher source term values is evident from the steep gradients observed in their surface temperature and surface pressure profiles. As a consequence, the temperature prediction of the present study appears better in the cavity closure region, while that of Hosangadi and Ahuja (2005) shows better agreement at the cavity leading edge. Thus, we emphasize again that the choice of model parameters poses a tradeoff, as we noticed in the global sensitivity analysis, in the prediction of pressure and/or temperature for cryogenic cases. Lastly, it is important to highlight that the inlet temperature gets closer to the critical temperature and the cavitation number decreases, as we proceed from case '283B' to '296B'. In comparison, the inlet velocity and consequently the Reynolds number assume values within the same order of magnitude for the depicted cases. Thus, under isothermal conditions, an increase in the cavity length is expected from case '283B' to '296B'. On the contrary, the surface pressure plots in Figure 22 clearly indicate a decrease in cavity length from case '283B' to '296B', despite the decrease in the freestream cavitation number. This fact clearly distills the significant impact of the thermal effect in cryogenic fluids, especially under working conditions that are close to the thermodynamic critical point. 290C Cavitation number 1.78 1.76 1.75 1.73 1.71 296B Cavitation number 1.78 1.74 1.70 1.67 1.63 1.62 Figure 23. Cavitation number (c = p py(T)/(0.5pU2)) based on the local vapor pressure Merkle et al. Model. Note the values of o = P. p,(TO)/(O.5pU2) for the cases '290C' and '296B' are 1.7 and 1.61, respectively. Our contention on the thermal effect is corroborated by the cavitation number (o = [p P(T)]/(0.5pU2) ) contours depicted in Figure 23. The freestream cavitation number (o.) of the case '296B' is smaller than the case '290C' (Table 5). However, the combination of evaporative cooling and its resultant impact over the vapor pressure causes a sharp increase in the effective cavitation number close to the cavitation zone. This increase is more substantial for the case '296B' and eventually leads to comparable levels of effective cavitation number between the two cases, as seen in Figure 23. (a) S089 069 050 S030 (b) 089 069 050 S030 (c) S089 0 69 0 50 0 30 Figure 24. Cavity shape indicated by liquid phase fraction for case '290C'. Arrowed lines denote streamlines (a) Merkle et al. Model isothermal assumption (b) Merkle et al. Model with thermal effects (c) Mushy IDM with thermal effects 63 Figure 24 reveals the phase fraction distribution and the streamlines for the computational cases of '290C'. The two models differ noticeably at the rear end of the cavitation zone. The cavity of the Mushy IDM consistently indicates lower condensation rates in contrast to the Merkle et al. Model, because of its longer length. We also note the gradual variation in density and the less extent of vapor phase in the cavities, which highlight the mushy/soggy nature of cavitation in cryogenic fluids. This is unlike our previous findings on regular fluids such as water (Senocak and Shyy 2004a, 2004b; Wu et al. 2003c). Under noncryogenic conditions, the flow structure in the cavitation vicinity is generally characterized by large streamline curvatures and formation of recirculation zones (Senocak and Shyy 2002, 2004a, 2004b; Wu et al. 2003c) (also seen in Figure 13). However, the weak intensity of cavitation in cryogenic fluids has a modest impact over the flow structure, as indicated by the streamline patterns in Figure 24. Figure 24(a) illustrates a 'special' solution to the case '290C' with the Merkle et al. Model assuming isothermal assumptions (energy equation not solved). Merkle et al. Model m f ( m 36E+02 2 37E+03 3 69E+03 3 95E+02 6 85E+03 2 50E+01 1 OOE+04 I 1 36E+00 Mushy IDM m ( m 1 536E+02 2 37E+03 369E+03 3 95E+02 685E+03 2 50E+01 1 OOE+04 1 36E+00 Figure 25 Evaporation (ih ) and condensation (i+ ) source term contours between the two cavitation models case '290C'. Refer to equations (3.8) and (3.18) for the formulations. Note that discounting the thermodynamic behavior yields a substantially large cavity size under similar conditions. Based on these observations, we emphasize that, from a modeling standpoint, the thermal effect is manifested via a combination of the cavitation model adaptations and the temperature dependence of physical properties. We bolster our argument on the condensation rates between the two models in Figure 25. Note that the evaporation and condensation contour plots demonstrate a mutually exclusive behavior. The condensation region of the Merkle et al. Model illustrates sharper gradients and is effective over a much larger region, as compared to the Mushy IDM. Following our assessment with liquid Nitrogen, we extend our focus to the cases with liquid Hydrogen, which has a density ratio of 0(30), unlike the 0(100) value for Nitrogen. We experience the need to recalibrate our cavitation models for this different fluid because of the discernible role played by the density terms (p,, p, and p,) in determining the volume transfer rates. Our numerical experimentation yields Cd,, = 0.82 and Cprod = 54.4 as appropriate values for the Merkle et al. Model, in case of liquid Hydrogen. Consistently, we choose / = 0.065 for the Mushy IDM in context of liquid Hydrogen. Our case selection for liquid Hydrogen, as seen from Table 6, follows similar trends as in case of liquid Nitrogen. The inlet velocity for all the cases with liquid Hydrogen is greater than 50 m/s, and it increases from the case '248C' to '255C' (66.4 m/s for '255C'). As a result, we upfront underscore these substantially higher values of inlet velocities that are employed for Hydrogen. Figure 26 depicts reasonable balance between the temperature and pressure predictions for the case '248C'; however, the agreement with the experimental data deteriorates equally for both the models as we proceed from the case '248C' to the case '255C'. Especially, the temperature is highly underpredicted as the inlet velocity increases between the two limiting cases.  Mushy IDM  Merkle et al. Model U Hord,1973;p 248C 0 05 1 x (inches from bend)  Mushy IDM  Merkle et al. Model U Hord,1973;p A Hosangadi &Ahuja, 2005 i 249D 05 0 05 1 x (inches from bend) 15  MushylDM  Merkle et al. Model U Hord,1973;p 255C 05 0 05 1 15 x (inches from bend)  Mushy DM  Merkle et al. Model U Hord,1973;T 248C 0 05 1 x (inches from bend) 15 2  Mushy DM  Merkle et al. Model U Hord,1973;T A Hosangadi &Ahuja,2005 249D 0 05 1 x (inches from bend) 15 2  Mushy IDM  Merkle et al. Model U Hord,1973;T S255C 0 05 1 15 2 x (inches from bend) Figure 26. Surface pressure and temperature for 2D hydrofoil for cases involving liquid Hydrogen. The results referenced as 'Mushy IDM' and 'Merkle et al. Model' are contributions of the present study. i I,,,,I,,,,I,,,,I I I II,, I I I I III, ,II, ,II, I The results of Hosangadi and Ahuja (2005), while depicting the consistent aspect of sharp gradients that we noticed earlier, also portray significant discrepancies between the surface temperature and the experimental data at the rear region of the cavity. This finding could be either attributed to inadequacies in the cavitation model parameters or the representation of temperaturedependent physical properties. However, the observed sensitivity of predictions to physical properties, the relatively better agreement in the case '248C', and the clear trend of disagreement with increasing velocities indicate the likelihood of the latter reason. Our extraction of physical properties from the NIST database (Lemmon et al. 2002) is based on models which assume thermodynamic equilibrium conditions. However, at such high fluid velocities, the timescales for thermodynamic equilibrium may be much larger than the overall flow timescales. This may introduce substantial discrepancies in the values of various physical properties and subsequently in the predictions; refer to Hosangadi and Ahuja (2005) for similar reporting. Of course, confirmation of the above possibility and development of rigorous nonequilibrium strategies is a challenging proposition for future research, and requires additional experimental insight. Nonetheless, we note the consistent performance of the Mushy IDM, especially in terms of the pressure prediction, over the chosen range of reference temperatures and velocities. We finally attempt to reassess our above calibrations (Cde, Cprod, and / ) for both the fluids and instill confidence into our predictive capabilities by performing computations over the axisymmetric 0.357inch ogive geometry. The ogive surface pressure and temperature plots for all the three cases are displayed in Figure 27. The shrinkage of the cavity length between cases '312D' and '322E', despite the decrease in 67 the reference cavitation number, unequivocally reemphasizes the influence of the thermosensible working conditions. 0 05 1 x in inches from bend 322E 0 02 04 06 08 1 x in inches from bend 349B 0 02 04 06 08 1 x in inches from bend  Mushy DM  Merkle et al. Model U Hord,1973; p 12 14  Mushy DM  Merkle et al. Model U Hord,1973; p 12 14 312D I F1  Mushy IDM  Merkle et al. Model U Hord,1973; T 0 05 1 x in inches from bend 322E  MushyDM  Merkle et al. Model U Hord,1973;T 4 0 05 1 x in inches from bend 349B 15  MushyDM  Merkle et al. Model U Hord,1973;T 0 05 1 x in inches from bend 15 Figure 27. Surface pressure and temperature for axisymmetric ogive for all the cases (Nitrogen and Hydrogen). The results referenced as 'Mushy IDM' and 'Merkle et al. Model' are contributions of the present study. i , I I I I I I I I I Furthermore, the discrepancy in our predictions for liquid Hydrogen, when subjected to higher velocities, is also evidenced by the plots for the case '349B'. Overall, the agreement with experimental data in all the cases is better in terms of the surface pressure than temperature. But, unlike the hydrofoil geometry, the temperature at the first (leading) probe point is overpredicted for all the ogive cases. As expected, this discrepancy is most in the liquid Hydrogen case ('349B'). Thus, the two geometries do not produce a consistent pattern of disagreement midst the surface temperature and the experimental measurements. This inconsistent behavior warrants further experimental and numerical probing from a standpoint of developing more precise cavitation models for cryogenic cavitation. Nonetheless, we observe that our analyses/calibrations are able to yield a justifiable range of results for the ogive geometry as well. CHAPTER 4 TIMEDEPENDENT COMPUTATIONS FOR FLOWS INVOLVING PHASE CHANGE Researchers have striven to develop efficient and accurate methodologies to simulate timedependent cavitating flows. The multiphase nature of the flow accompanied by complex flow physics yields a system of tightlycoupled governing equations. Furthermore, interfacial dynamics, compressibility effects in mixture region, and turbulence entail deployment of numerical models to represent these physical phenomena. The formulation of these models substantially impacts the solution procedure. As a result, evolving efficient algorithms for unsteady cavitating flows is certainly a nontrivial task. Formulation of implicit procedures for pressurebased methods is impeded mainly by the strong linkage between the flow variables such as velocity and pressure. Furthermore, the dynamics of the variable density field in cavitating flows imposes supplementary equations on the existing system, which add to the computational challenges. As a result, iterative algorithms such as SIMPLE, SIMPLER, and SIMPLEC (Versteeg and Malalasekera 1995), which are commonly employed for a wide variety of problems, may be computationally expensive for solution of cavitating flows with pronounced unsteady behavior. Senocak and Shyy (2002, 2004b) circumvented this difficulty by incorporating the Pressure Implicit with Splitting of Operators (PISO) algorithm suitably with the model equation of cavitation dynamics and the sound propagation model for the mixture region. This endeavor resulted in a noniterative methodology for timedependent computation of cavitating flow through a series of predictorcorrector steps. Senocak and Shyy (2002, 2004b) demonstrated the merits of this efficient approach with a series of cavitating flow computations on geometries such as ConvergentDivergent Nozzle and hemispherical solids. However, all the past efforts were in the context of isothermal cavitation sans the energy implications and real fluid properties. From standpoint of formulating a non iterative algorithm for cryogenic cases, the inclusion of the nonlinear energy equation and temperature dependence of physical properties into the existing methodology is expected to pose multitude of adversities to the solution efficiency and accuracy. Direct development of an algorithm for cryogenic cavitation may be a presumptuous preliminary approach. As a result, we initiate a multiphase noniterative procedure on a simpler test problem with similar nature of challenges, and, with an objective of monitoring the accuracy and stability of the computations. In this chapter, we elucidate the newly initiated algorithm and illustrate its results on a chosen test case. The truly unsteady problem of Gallium fusion (with natural convection effects) is adopted for the purpose of validation. Discussion on the grid sensitivity, accuracy of results, and the stability criterion is provided. The above primary objective is complimented by a reducedorder description of the gallium fusion problem by Proper Orthogonal Decomposition (POD). The flowfield in the problem is characterized by a solidliquid front movement with a continuous change in the overall flow length scale. The ability of POD, to accommodate these varying flow scales, is mainly emphasized. The insights gained from this study on gallium fusion are later extended to improving the algorithm for cryogenic cavitation. The algorithm in context of cryogenic cavitation is described to point out the changes implemented in order to address the temperature effects. POD is employed to probe the timedependent data, and offer it a succinct representation. 4.1 Gallium Fusion Experimental studies (Gau and Viskanta 1986) have investigated the physics of Gallium fusion due to its low fusion temperature and ease of handling. The availability of 2D experimental data has motivated numerical studies (Lacroix 1989, Lacroix and Voller 1990, Shyy et al. 1995) to adopt it as a test case. T=0 Y=l 7 0 liquid solid T=1 T=O Y _X aT = 0 X= 1 X y Figure 28. Schematic of the 2D Gallium square geometry with the Boundary Conditions (Shyy et al. 1998) A schematic of the square geometry along with the boundary conditions can be viewed in Figure 28. 4.1.1 Governing Equations The singlefluid modeling approach is adopted by employing a liquid phase fraction f The density is assumed constant through the Buossinesq approximation. The governing equations for the problem are as follows. au = 0 (4.1) O(put) O(puuox) ap a au 1 f2 + + (P )C(f )u +gpf(T T)2 (4.2) at axi ax, Qxj cx f3 +q a(ph) a a aT + (puh)= (k ) (4.3) at )x ax ax Here, / represents the coefficient of thermal expansion for the fluid. The term 1 f2 C( )u, in the momentum equations is the Darcy source term (Shyy et al. 1998). f +q Through their functional dependence on f they are modeled to retard the velocity to insignificant values in solid region. The constants C and q, in compliance, are tuned to yield a negative source term, which is at least seven orders of magnitude higher than other terms, in the solid region (f= 0). The enthalpy in equation (4.3) is expressed as: h= CT + jL (4.4) The phase fraction f is modeled computationally by the hbased method (Shyy et al. 1998). The enthalpy at any given iteration is computed as shown below followed by an update off(explained in a later section). hk =CpTk +fkL (4.5) Note that f is always bounded by 0 and 1. Thus, equation (4.5) can be recast into a temperature equation as shown. a(pCPT) a a aT a(pLf) a CT) (pCuT)= (k )[O+ (pLujf)] (4.6) at 8x) 8x^ 8x) at axi The iterative update of the liquid phase fraction, as mentioned above, renders the energy/temperature equation nonlinear. Furthermore, the buoyancy term in equation (4.2) leads to a pressurevelocitytemperature coupling. The nonlinearity in the energy equation and the strong coupling between various flow variables justify the use of the gallium fusion problem as a test case preceding the cryogenic problem. The physical properties such as viscosity, latent heat, thermal conductivity, and specific heat are, however, assumed to be constant during the fusion process. 4.1.2 Numerical Algorithm The previous computations on similar multiphase problems employed an iterative solution strategy (Chuan et al. 1991, Khodadadi and Zhang 2001). Conversely, the PISO method (Issa 1985, Thakur et al. 2002, Thakur and Wright 2004), which forms the backbone of the current algorithm, essentially is a series of predictor corrector steps to yield a noniterative solution of flow equations through operator splitting procedure. The present methodology, however, closely follows a slightly modified version of PISO designed for buoyancy driven singlephase flows (Oliveira and Issa 2001). Kim et al. (2000) proposed a series of steps to accelerate the simple heat conduction equation with a solidliquid phase boundary. The following methodology attempts to blend the merits of the modified PISO (Oliveira and Issa 2001) with the propositions of Kim et al. (2000), to design an efficient and accurate algorithm for phase change problems. The sequence of calculations for the algorithm, based on the above governing equations, is expounded below. The strongly implicit form of the discretized governing equations with a finite volume formulation, at any node P, is as follows: A (pu, ) = 0 (A)"u =_ A1"u) Ap1+ +H"up +Spu(un+1) (APv +1 = A vb'1 A pn+l + H;v + S (vn1) (4.7) +Bp (T n+) ) TTn+1 TT;II T n n + M (fn (AP)T + = A '+HT +Hfp +M(f ) Here, A represents the terms at the current time level, H represents the terms at older time step, M represents the terms involving the phase fraction f B represents the buoyancy term, and S represents the Darcy source term. The Darcy source terms can be absorbed into the left hand side term to yield to following equations. A (pul n+') 0 (G;P1 )p' = C AnInln A "n+l + Hiu" (4.8) (Gpv n)1 = ZAb1 A, pn1 +Hvp + B (Tn+1) STn+l T T +I T n n + + ( n (AP)T' = A b +HP +H f/ +M(f"1) Here, Gp" = Ap +Sp" and so on. The sequence of steps to solve equation (4.8) non iteratively is elaborated next. (a) Momentum predictor (GPu)uP = Au au" + H u (4.9) (Gp)v = ZA v A avp" + H + BPv (T+ ) These velocities u* and v* are obtained by using the pressure value at the previous time step and hence are not divergence free. (b) First pressure corrector A, (pu,) = 0 p (4.10) A A, p')= A,(pu) ) GP Here, p'=p*p is calculated. The pressure field p* can now be employed for correcting the divergence error in velocity. Note that the pressure correction at this stage is limited by several approximations and is not adequate to produce an accurate velocity field. (c) First momentum corrector (GPu)u**= _Au Ap +Hu"p (4.11) (GP )vP *= 'Avn A p* + Hv + B(T") The velocities are corrected to yield u** andv** explicitly using the intermediate pressure field obtained in the previous step. (d) First temperature corrector (AP )TP = Ab + +HTP" +Hf fP +M(f*) (4.12) The above temperature equation utilizes the latest value off However, in view of a non iterative strategy, solution of temperature, merely by the above equation, may not be sufficient for rapid convergence due to the delayed update off As a remedy, a series of explicit steps are implemented following the above equation to significantly improve the prediction of temperature and, subsequently, f These steps are adapted from the algorithm proposed by Kim et al. (2000) for pure conduction equation. They are enlisted below. (i) Updatefby hbased method as shown. hk+1 k+1 k fL f = 0 if hk+l < f~+l =1 if hk > h, (4.13) f, k P elsewhere h,hk where, h = CpT,h, = CpT, +L Note that at the outset T7' = Tp. (ii) If fpl fp correct temperature using the explicit form of equation (4.12) with updated coefficients. This step yields the new temperature field T+2 (iii) Compute residual of the temperature equation as shown further. O = (AP)T A, +HT +Hf f/ +M(fk+l) (4.14) The above residual equation can be further employed to improve temperature prediction by NewtonRaphson approach as follows (Kim et al. 2000): aTk+2 P where, the Jacobian is approximated as 0+ 1 Af . P P In summary, steps (i) (iii) are performed at least 1012 times following the equation (4.12). Let the temperature and phase fraction field at the end of this stage be denoted as T* andf*. (e) Second pressure corrector A (pul ) = 0 SAP") = ( Bp(T*)Bp(T")+ G"p Gp" (4.16) SAbu** Abu* [(SP)* (SP)]** }) This step is a powerful characteristic of the PISO algorithm. The 'differential Darcy term', (Sp')* (Sp)up**, is an addition to the terms already suggested by Issa (1985) and Oliveira and Issa (2001). This term arises out of the update of G" after the first temperature corrector. Note that the coefficient G"' unlike singlephase flow cases, also includes the implicit coefficient of the Darcy term (refer to the discussion in step (a)). The Darcy term undergoes a sudden change over several orders of magnitude between the liquid and the solid region. Since the first temperature corrector tends to change the phase fraction distribution, it is imperative to update the coefficient Gp post that step. The second pressure corrector, which is obtained via subtracting equation (4.11) from (4.17), is thus able to account for phase front movement because of the 'differential Darcy term'. In summary, the second pressure corrector attempts to couple the nonlinear terms in the momentum equations, the temperature field, and the phase front movement to the pressure field. (f) Second momentum corrector (Gpu*)u* = C Auu u p Ap** + Hu *** ** ** (4.17) (Gpv )v = Avb Ap + Hv"+Bv (T) The velocity field at this stage is derived from the more accurate prediction of pressure p** obtained from the previous step. (g) Second temperature corrector (A)'T* T** +HT+H f, +AM*(f**) (4.18) The prime purpose of this equation is to correct the temperature field with updated velocities. This step is followed by an update off, similar to that in equation (4.13). It was reported that repeating the corrector steps (steps (b) (g)) for an extra time improves coupling between various equations (Thakur and Wright 2004). However, due to the large magnitude of nonlinearity expected for the fusion problem, at least 3 more repetitions of steps (b)(g) are recommended (totally 8 corrector steps). The above modification to the PISO, similar to the original algorithm, incurs a splitting error due to the operator splitting approach. However, the order of this splitting error is higher than the formal order of temporal accuracy (Issa 1985), thus justifying the computational accuracy. 4.1.3 Results The multiphase algorithm is implemented in conjunction with a pressurebased solver having multiblock capability (Thakur et al. 2002). The diffusive terms are handled by central differencing, while the first order upwind scheme is employed for the convective terms. Computations are mainly performed on a 2D square domain of size D=1. The range of values for St and Ra in the computations is [1 0.042] and [104 2.2x106] respectively. Sample results obtained from current algorithm are compared to those obtained by Shyy et al. (1995) in the following discussion. It is important to note that results obtained by Shyy et al. (1995) have greater fidelity than those published in similar studies, due to the use of a finer grid. 4.1.3.1 Accuracy and grid dependence Figure 29. 2D interface location at various instants for St = 0.042, Ra = 2.2 x 105 and Pr = 0.0208. White circles represent interface locations obtained by Shyy et al. (1995) on a 41x41 grid at time instants at t = 56.7s, 141.8s, & 227s respectively. Figure 29 illustrates the 2D interface locations obtained with the current algorithm, with three different grid resolutions, for the parameters shown. Note that some results by Shyy et al. (1995) are at a slightly earlier time instant. This contributes to the modest difference observed in the interface location, especially during the early stages, when the interface velocity is significantly higher. Considering this fact, the present results show reasonable agreement with the earlier computation. Furthermore, the timedependent movement of the interface is consistent on all the three grids, although there is a qualitative impact of the resolution on the interface profile. 0t5 0 0 x 0 x0    21x2l grid (a) 21x21 gdd 41x41 gdd S 81x81gdd ; ; S 025 05 075 X (b) 1 (c) x x Figure 30. Grid sensitivity for the St = 0.042, Ra = 2.2 x 105 and Pr = 0.0208, 2D case (a) Centerline vertical velocity profiles at t = 227s (b) Flow structure in the upper left domain at t = 57s; 41x41 grid (c) Flow structure in the upperleft domain at t= 57s; 81x81 grid. Table 7. Location of the primary vortex for the St = 0.042, Ra = 2.2 x 105 and Pr = 0.0208 case Grid size xlocation at t ylocation at t xlocation at t ylocation at t = 57s = 57s = 227s = 227s 21x21 0.130 0.760 0.350 0.674 41x41 0.135 0.755 0.365 0.661 81x81 0.146 0.710 0.377 0.658 The influence of the grid quality on the flow structure can be assessed from the centerline velocity profiles and the locations of the center of the primary vortex depicted in Figure 30(a) and Table 7 respectively. The even number of nodes (cell centers) in each grid creates a slight offset between the center gridline and the geometric centerline. This offset varies inversely with the grid quality, and is expected to contribute to the reasonable discrepancy in the centerline velocity profiles, shown in Figure 30(a). In comparison, the center of the primary vortex, from the findings reported in Table 7, indicates a gradual, downward shift with grid refinement, particularly during the initial stages of fusion. In fact, this movement, although moderate, is conspicuous at t = 57s, when an 81point grid is used instead of a 41point grid. This observation has a physical relevance, as noticed from Figure 30(b) and Figure 30(c). The upperleft part of the domain develops a secondary vortex after some initial time lapse. This vortex is expected to induce a downward motion of the primary vortex structure. As clearly seen, the 41x41 grid, unlike the finer grid, is unable to capture this smallscale circulation, which justify the data in Table 7. The overall flow convection pattern and, consequently, the interface movement are, however, weakly affected by the secondary vortex structure. Furthermore, the solution demonstrates a fairly consistent improvement with the grid quality, in addition to a restrained griddependence. These findings are encouraging from the standpoint of stability restrictions, which are elucidated in the following section. Although the multiphase algorithm is elucidated/illustrated in context of 2D calculations, its generality for 3D computations is examined by extending the case in Figure 29 to square box geometry. The fusion process is initiated by two adjacent, heated walls to yield spanwise flow variations. 