COMPUTATIONAL METHODOLOGY FOR THE SIMULATION OF TURBULENT CAVITATING FLOWS By NAN ENOCAK A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2002
Copyright 2002 by nan enocak
To my parents and my brother
ACKNOWLEDGMENTS I would like to express my sincere thanks and appreciation to my advisor Professor Wei Shyy for his guidance and support throughout this research. His tolerant and generous attitude helped me to accomplish this task. I also would like to thank Professors Louis N. Cattafesta, James F. Klausner, Renwei Mei, and Mark Sheplak for their valuable comments and for serving on my committee. My thanks also go to all the individuals of the Computational Thermo-Fluids group, whom I have the privilege to work with. In particular, Marianne Franois has been a close friend over the past few years. My parents have always set their sonsâ€™ education as their top priority. I have benefited from this attitude throughout my education. I am thankful beyond words. My brother, Kvan enocak, has been a great support during my study. Together we had memorable and joyful times. Finally, I would like to thank National Science Foundation (NSF) and Office of Naval Research (ONR) for the financial support. iv
TABLE OF CONTENTS page ACKNOWLEDGMENTS.................................................................................................iv LIST OF TABLES............................................................................................................vii LIST OF FIGURES.........................................................................................................viii LIST OF SYMBOLS........................................................................................................xii ABSTRACT.......................................................................................................................xv CHAPTER 1 INTRODUCTION.......................................................................................................1 1.1 Research Motivation and Objectives................................................................1 1.2 Contributions of the Present Study...................................................................3 2 LITERATURE REVIEW............................................................................................4 2.1 Fundamentals of Cavitation..............................................................................4 2.1.1 Types of Cavitation....................................................................................4 2.1.2 Effects of Cavitation..................................................................................6 2.1.3 Thermodynamics of Cavitation..................................................................7 2.1.4 Interfacial Conditions and Bubble Dynamics............................................9 2.1.5 Cavitation Inception and Cavitation Parameter.......................................14 2.1.6 Kinetic Theory of Mass Transfer.............................................................15 2.2 Experimental Studies......................................................................................17 2.3 Numerical Studies...........................................................................................23 3 THEORETICAL FORMULATION..........................................................................27 3.1 Governing Equations......................................................................................27 3.2 Turbulence Modeling......................................................................................30 3.2.1 Two Equation (k-) Turbulence Model...................................................31 3.2.2 Wall Functions.........................................................................................32 3.2.3 Implications of the Near-wall Turbulence Treatment..............................33 3.3 Cavitation Modeling.......................................................................................38 3.3.1 Transport Equation-Based Empirical Cavitation Models........................39 v
3.3.2 Derivation of a Transport Equation-Based Interfacial Dynamics Cavitation Model.....................................................................................41 3.4 Numerical Methodology.................................................................................47 3.4.1 Pressure-Based Algorithm for Steady-State Computations of Cavitating Flows......................................................................................50 3.4.2 Pressure-Based Algorithm for Time-Dependent Computations of Cavitating Flows......................................................................................56 3.4.3 The Issue of Speed of Sound Definition..................................................60 3.4.4 Boundary Conditions...............................................................................62 4 RESULTS OF STEADY-STATE COMPUTATIONS.............................................64 4.1 Validation of the Pressure-Based Algorithm for Cavitating Flows................64 4.1.1 Flow over a Hemispherical Object..........................................................65 4.1.2 Flow over a Blunt Object.........................................................................72 4.2 Evaluations of Transport Equation-Based Cavitation Models........................77 4.2.1 Flow over a Hemispherical Object..........................................................78 4.2.2 Flow over a NACA66MOD Hydrofoil....................................................82 4.2.3 Convergent-Divergent Nozzle (Stable Sheet Cavitation)........................88 5 RESULTS OF TIME-DEPENDENT COMPUTATIONS........................................96 5.1 Single-Phase Flow over a Cylinder.................................................................96 5.2 Unsteady Cavitating Flows in Convergent-Divergent Nozzles......................99 5.2.1 Implications of the Speed of Sound Definition......................................101 5.2.2 Cavitation Instability (Auto-Oscillations)..............................................105 5.2.3 Two-Phase Flow Structure.....................................................................109 6 CONCLUSIONS AND FUTURE DIRECTIONS..................................................114 6.1 Conclusions...................................................................................................114 6.2 Future Directions..........................................................................................117 LIST OF REFERENCES................................................................................................118 BIOGRAPHICAL SKETCH..........................................................................................125 vi
LIST OF TABLES Table page 2-1 Physical properties of water at standard atmospheric pressure..................................9 2-2 An overview of the past experimental studies on cavitating flows..........................22 2-3 An overview of the past numerical studies on cavitating flows...............................26 3-1 The empirical constants used in the variants of kturbulence model.....................32 3-2 Description of grid properties...................................................................................34 5-1 Influence of top and bottom boundary conditions on the Strouhal number.............97 vii
LIST OF FIGURES Figure page 2-1 Visualization of different types of cavitation.............................................................6 2-2 Surface erosion on ship propellers due to cavitation..................................................7 2-3 Pressure vs. specific volume diagram explaining the vaporization of a liquid..........8 2-4 Schematic of a liquid-vapor interface.......................................................................10 2-5 Schematic of a bubble in an infinite domain of liquid.............................................11 2-6 Schematic of the mass transfer through a cavity......................................................15 3-1 Grid distribution along j direction............................................................................35 3-2 The effect of grid distribution on wall shear stress..................................................35 3-3 y + profiles along the flat plate. =0........................................................................36 3-4 Effect of different y + definitions on wall shear stress. =0....................................37 3-6 Effect of different y + definitions on wall shear stress. =5....................................38 3-7 Representation of a vaporous cavity in homogeneous flow theory..........................42 4-1 The computational domain and imposed boundary conditions................................66 4-2 Grid refinement study for hemispherical object with =0.40..................................66 4-3 Improvement in convergence level due to density upwinding. Hemispherical object at =0.40........................................................................................................67 4-4 Effect of convection schemes on pressure coefficient and density distribution for the hemispherical object at =0.40.................................................................................69 4-5 Sensitivity of modeling parameters for the hemispherical object at =0.40............69 viii
4-6 Comparison of pressure coefficient distributions for hemispherical object under noncavitating and cavitating conditions....................................................................70 4-7 Vorticity generation at the closure region due to baroclinic torque.........................71 4-8 Dynamics of the phase change.................................................................................72 4-9 The computational domain and imposed boundary conditions................................73 4-10 Impact of turbulence modeling on predictions of pressure coefficient distribution for noncavitating conditions................................................................74 4-11 Effect of turbulence model on separation zones for noncavitating conditions.......75 4-12 Impact of turbulence modeling on predictions of pressure coefficient and density distribution.................................................................................................75 4-13 Comparison of u-velocity along the surface obtained from different turbulence models....................................................................................................................76 4-14 Comparison of turbulence quantities obtained from different turbulence models..77 4-15 Comparison of surface pressure distribution obtained from empirical cavitation models...................................................................................................79 4-16 Comparison of surface density distribution obtained from empirical cavitation models...................................................................................................80 4-17 Comparison of surface pressure and density distribution obtained from empirical Model-1 and the new interfacial dynamics model................................................81 4-18 Detailed view of the cavitating region. On the left is the density distribution; on the right is the spanwise vorticity distribution...................................................81 4-19 Visual comparison of the grid distributions in the cavitating region.....................83 4-20 Effect of grid distribution on predictions. Model-2 is used....................................83 4-21 Distribution of the y + values of the first grid points away from the wall along the hydrofoil surface..............................................................................................84 4-22 Comparison of surface pressure distribution. The cavitation number is 0.91.......86 4-23 Comparison of surface pressure distribution. The cavitation number is 0.84........87 4-24 Detailed flow visualization. Cross sectional view shows pressure distribution.....87 ix
4-25 Comparison of the present results obtained from Model-1 with the results of published studies in literature that also utilize Model-1.........................................88 4-26 The computational domain and the imposed boundary conditions for the convergent-divergent nozzle..................................................................................89 4-27 The impact of cavitation on the structure of internal flow through the nozzle......90 4-28 The effective terms of the vorticity transport equation for cavitating conditions..92 4-29 The velocity field around sheet cavitation..............................................................93 4-30 Comparison of the cavity boundary and the velocity profiles within the cavity....94 4-31 Vapor volume fraction profiles within the cavity...................................................95 5-1 Computational domain and imposed boundary conditions......................................97 5-2 Periodic oscillation of the cross-stream (v) component of velocity.........................98 5-3 Spanwise vorticity distribution. Snapshot of the vortex shedding and the von Karman vortex street is shown.................................................................................98 5-4 Nozzle-1, computational domain and the imposed boundary conditions.................99 5-5 Nozzle-2, computational domain and the imposed boundary conditions...............100 5-6 The impact of the speed of sound definition on time-dependent behavior. The new interfacial dynamics cavitation model is used. Geometry is Nozzle-1....102 5-7 Comparison of the time-dependent behavior of different cavitation models. Geometry is Nozzle-2............................................................................................103 5-8 The impact of the speed of sound definition on time-dependent behavior. Geometry is Nozzle-2. The cavitation number is 0.67..........................................104 5-9 Cavity shapes and the recirculation zone during the instability cycle. Geometry is Nozzle-1. The new interfacial dynamics cavitation model is used........................108 5-10 The periodic oscillation of the upstream pressure in Nozzle-2. The cavitation number is 1.98 based on time-averaged quantities...............................................109 5-11 Time averaged velocity profiles (u/U ref ) within the cavity formed in Nozzle-1..110 5-12 Vapor volume fraction profiles within the cavity formed in Nozzle-1................111 x
5-13 Time-averaged velocity profiles (u/U ref ) within the cavity formed in Nozzle-2..112 5-14 Time-averaged vapor volume fraction profiles within the cavity formed in Nozzle-2...............................................................................................................113 xi
LIST OF SYMBOLS b local characteristic speed c speed of sound C arbitrary O(1) constant C 1 , C 2 turbulence model constants C dest empirical constant in the evaporation term C turbulence model constant C P pressure coefficient C prod empirical constant in the condensation term E log-law constant f flux at a computational cell face F flow rate at a cell face h latent heat of vaporization j flux of molecules k thermal conductivity in Chapter 2, turbulent kinetic energy m evaporation rate m condensation rate M molecular weight L ch characteristic length P pressure, turbulence production in Eq. (3-15) q source term r radial coordinate R radius of curvature u velocity vector u i velocity components in Cartesian coordinates u + nondimensional velocity u * friction velocity xii
U magnitude of the horizontal component of velocity t, t time, mean flow time scale T temperature, mean flow time scale in Chapter 5 x i Cartesian coordinates y + nondimensional normal distance from the wall y P normal distance of the first grid point away from the wall volume fraction ij Kronecker delta function turbulent dissipation rate v pressure increment in m surface tension diffusion coefficient log-law constant laminar viscosity eff effective viscosity t turbulent viscosity , , curvilinear coordinates generalized dependent variable ij shear stress Rij Reynolds stress wall wall shear stress t kinematic viscosity mixture density cavitation parameter k , turbulence model constants Subscripts, Superscripts B bubble cf computational cell face xiii
I interface L, l liquid phase V,v vapor phase m mixture phase n normal direction s tangential direction t tangential direction x component in x coordinate direction y component in y coordinate direction z component in z coordinate direction free stream * predicted value xiv
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 METHODOLOGY FOR THE SIMULATION OF TURBULENT CAVITATING FLOWS By nan enocak August 2002 Chair: Wei Shyy Department: Aerospace Engineering, Mechanics and Engineering Science A comprehensive computational methodology is developed for turbulent cavitating flows. The ensemble-averaged Navier-Stokes equations, along with a volume fraction transport equation, are employed. The flow field is computed in both phases with the vapor pressure recovered inside the cavity via a mass transfer model. The two-equation kmodel is adopted for turbulence closure. Based on the mean flow variables, both steady-state and transient computations are conducted. To ensure stable and convergent computations, unified incompressible-compressible pressure-based algorithms for steady and time dependent computations are presented, with the large density ratio at the phase boundary handled via pressure-density coupling schemes. While no temperature effect is considered, the resulting pressure-correction equation shares common features with that of high-speed flows, exhibiting a convective-diffusive nature instead of only a diffusive type. To address the empiricism in existing cavitation models, an analysis of the mass and momentum conservation at a phase change boundary is conducted and a new transport xv
equation-based interfacial dynamics cavitation model is developed. The new model eliminates the empirical parameters in existing models by introducing interfacial velocity terms. The computational capability has been applied to axisymmetric projectiles and a planar hydrofoil for steady-state simulation, and to convergent-divergent nozzles for both steady-state and time-dependent simulation. The steady-state results indicate that comparable pressure distribution predictions are yielded by all the cavitation models considered. The vorticity dynamics in the closure region is detailed and it is shown that the baroclinic vorticity generation is important for the turbulence production in the cavity closure region. To account for the two-phase mixture effect, the implications of the speed of sound definition in cavitating flows are investigated. The interplay of cavitation number, pressure gradient and cavitation model has been probed to assess the time-dependency of the cavitation dynamics. The computational methodology has demonstrated capabilities to predict both qualitative and quantitative features of the turbulent cavitating flow. xvi
CHAPTER 1 INTRODUCTION 1.1 Research Motivation and Objectives In liquid flows, cavitation generally occurs if the pressure in certain locations drops below the vapor pressure and consequently the negative pressures are relieved by the formation of gas-filled or gasand vapor-filled cavities (Batchelor 1967). Cavitation can be observed in a wide variety of propulsion and power systems like pumps, nozzles, injectors, marine propellers, hydrofoils and underwater bodies (Knapp et al. 1970). Cavitating flows in most engineering systems are turbulent. Dynamics of the interface formed involves complex interactions between vapor and liquid phases. These interactions are not well understood in the closure region of the cavity where a distinct interface may not exist and the flow is unsteady. Undesirable features of cavitation are structural damage, noise and power loss. On the other hand, drag reduction can be observed on bodies surrounded fully or partially with a natural or gas-ventilated cavity (Lecoffre 1999). One well-known earlier study on cavitation dates back to the research of Rayleigh (1917). The collapse of an empty cavity in a large mass of liquid was studied. Since then, a considerable amount of theoretical and experimental research has been conducted. Cavitation has been recognized for its damage to hydrodynamic structures and accordingly treated as a condition to be avoided. To eliminate these negative effects, both experimental and computational studies have been undertaken. To date, the efforts to 1
2 predict and model cavitation have been driven mainly by the turbomachinery industry. A computational model that predicts the location and extent of the cavity is often sufficient to be incorporated into design iterations of turbomachinery components. However, with improved analysis and design capability, more quantitative information, such as pressure distribution and time-dependent nature of the fluid flow are of growing interest. Apart from the damage associated with cavitation, there are applications in which one can benefit from cavitation. High-speed hydrofoil boats, supercavitating propellers and supercavitating vehicles are some examples. By generating supercavitation around hydrodynamic vehicles skin-friction drag can be decreased substantially. The gaseous cavity surrounding the vehicle can also provide a mechanism for improved control and guidance strategies. To achieve these goals, accurate and detailed information of the flow field is required. To date, a comprehensive computational model has not been established that can provide the detailed description of the turbulent cavitating flows. In cavitating flow modeling, regions of vapor phase are not known a priori. This implies that the major goal of a cavitation model should be to predict the onset, growth, departure break-up and collapse of bubbles in cavitating flows. Existing cavitation models compute the overall behavior of cavitating flows (such as pressure distribution on solid surfaces and cavity shape). There is no comprehensive model in the literature that can simulate various types of cavitation and provide a detailed description of the flow field.
3 The major objectives of this research are as follows: Develop a computational methodology that can provide a detailed description of turbulent cavitating flows and be able to simulate steady and unsteady cavitation. Investigate the physics of cavitating flows to provide improved understanding of the phenomenon based on the developed computational methodology. 1.2 Contributions of the Present Study The major contributions of the present study are summarized below. Investigation of the near wall turbulence treatment with the original kmodel and wall functions (He et al. 2000). A review of the computational studies on cavitation (Wang et al. 2001). Development of a unified incompressible-compressible pressure-based method for turbulent cavitating flows that can handle large density jump at the phase boundary (Senocak and Shyy 2002). Development of a transport equation-based interfacial dynamics cavitation model (Section 3.3.2) and assessment of the transport equation-based cavitation models in the context of the unified pressure-based algorithm. Development of a unified pressure-based operator splitting method for time-dependent computations of turbulent cavitating flows (Section 3.4.2). Identification of the importance of speed of sound definition in unsteady computations of cavitating flows (Section 5.2.1). An approximation to the fundamental definition of speed of sound (Section 3.4.2). Identification of the coupling of turbulence production and baroclinic vorticity generation downstream of cavitation (Section 4.2.3). Identification of the onset of cavitation instability (Section 5.2.2).
CHAPTER 2 LITERATURE REVIEW A fundamental understanding of cavitation is presented before summarizing the published experimental and numerical studies of cavitation. The literature review is centered on numerical and experimental studies that support the objectives of the present study. Among the numerous experimental studies, those that address the complex flow structure have been reviewed. 2.1 Fundamentals of Cavitation 2.1.1 Types of Cavitation Different types of cavitation can be observed depending on the flow conditions and geometry. Each type of cavitation has distinct characteristics. Five major types of cavitation have been described in the literature. In the following, each of them is described briefly. Traveling cavitation is a type of cavitation in which individual transient cavities or bubbles form in the liquid and move with it as they expand or shrink during their life cycles (Knapp et al. 1970). It is typically observed on hydrofoils at small angle of attack. The geometries of the bubbles formed are highly dependent on the amount of nuclei present in the incoming flow (Lecoffre 1999). To the naked eye traveling cavitation may appear as sheet cavitation. Individual transient spherical bubbles can be distinguished by snapshot photographs. Traveling cavitation on a hydrofoil is shown in Figure 2-1a. 4
5 Cloud cavitation is caused by vorticity shed into the flow field. It causes strong vibration, noise and erosion (Knapp et al. 1970). The shedding of cloud cavitation is periodic and the re-entrant jet is the basic mechanism to generate cloud cavitation (Kawanami et al. 1997). The process is shown in Figure 2-1b. Sheet cavitation is also known as fixed, attached, cavity or pocket cavitation. Sheet cavitation is stable in a quasi-steady sense. The phenomenon is shown in Figure 2-1c. The interface between the liquid and the vapor can be smooth and transparent or it can have the shape of a highly turbulent boiling surface (Knapp et al. 1970). The liquid vapor interface becomes wavy and breaks down in the closure region of the cavity. Downstream flow, which contains large scale eddies, is dominated by bubble clusters (Gopalan and Katz 2000). Supercavitation is when the sheet cavity grows in such a way to envelope the whole solid body (Knapp et al. 1970). Ventilation can be used to create or to enhance a supercavity. Supercavitation is desirable to achieve viscous drag reduction on underwater vehicles operating at high speeds. High speeds in excess of 1500 m/s are reported for supersonic operation of underwater projectiles (Kirschner 2001). Other than viscous drag reduction, lift force acting on hydrofoils can be increased by creating a supercavity on the upper surface at a constant pressure. A typical supercavitating hydrofoil is shown in Figure 2-1d. The vortex cavitation occurring on the tips of rotating blades is known as tip vortex cavitation (Figure 2-1e). Cavities form in the cores of vortices in regions of high shear. This type of cavitation is not restricted to rotating blades. It can also occur in the separation zones of bluff bodies (Knapp et al. 1970).
6 (a) (b) (c) (d) (e) Figure 2-1. Visualization of 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 the permission of EDP Sciences). 2.1.2 Effects of Cavitation Cavitation has both desirable and undesirable features. Undesirable effects are usually associated with surface erosion, excessive noise generation and hydrodynamic losses.
7 The collapse of cavitation bubbles is violent and leads to surface erosion and noise generation. Figure 2-2 shows the substantial damage on a ship propeller because of operation under cavitating conditions. Formation of cavities interrupts the interaction of the solid boundaries with the liquid phase. In internal flow devices, cavitation leads to a reduction in flow rate. Usually in pumps, considerable head losses occur because of cavitation. For instance, in pumps the occurrence of cavitation decreases the force that can be applied on the liquid (Knapp et al. 1970). This changes the angular deflection of the liquid, and as a result the pump operates at off-design conditions. Desirable features of cavitation lead to useful applications like micro-bubble generation, cleaning of objects, catalysts of chemical reactions and viscous drag reduction (Lecoffre 1999). Figure 2-2. Surface erosion on ship propellers due to cavitation. (Photograph by Inanc Senocak) 2.1.3 Thermodynamics of Cavitation Vaporization of a liquid can be through the addition of heat (in the case of boiling) or through the depressurization of the liquid (in the case of cavitation). Boiling and cavitation differs in certain aspects. Growth and collapse of bubbles are observed in
8 cavitation, whereas continuous growth of bubbles is typical in boiling. These two modes of vaporization differ in the pressure vs. specific volume diagram also. As seen from , boiling follows an isotherm and the phase change occurs because of heat addition at constant pressure. On the other hand a shift in isotherms is observed during cavitation. Figure 2-3 Figure 2-3. Pressure vs. specific volume diagram explaining the vaporization of a liquid (adapted from Lecoffre 1999). vaporliquid Critical point boiling cavitation Pv isotherms vaporliquid Critical point boiling cavitation Pv isotherms Although thermal effects are usually neglected for liquids like water at room temperature, these effects may become important for cryogenic fluids. Despande et al. (1997) have computed the thermal boundary on the cavity surface. A temperature depression zone has been observed that resulted in the shortening of the cavity length. The physical properties of water widely encountered in cavitation research are shown in . Table 2-1
9 Table 2-1. Physical properties of water at standard atmospheric pressure. (adapted from Incropera and Dewitt 1996, Reynolds 1979) Saturated Water at 300 K L 997.0 Density (kg/m 3 ) v 0.0256 L 8.55x10 -4 Absolute viscosity (N.s/m 2 ) v 9.09x10 -6 k L 0.613 Thermal conductivity (W/m.K) k v 0.0195 Vapor pressure (kPa) P v 3.531 Surface Tension (N/m) 7.17x10 -2 h L 111.7 Enthalpy (kj/kg) h v 2550.1 2.1.4 Interfacial Conditions and Bubble Dynamics 18.104.22.168 Governing equations at a phase-change interface Consider an interface separating a liquid and its vapor during an evaporation process as shown in Figure 2-4. At the interface, the system must satisfy the conservation of mass, momentum and energy (Carey 1992). In the following formulation the velocities are with respect to a stationary observer. The normal velocity of the liquid moving toward the interface is represented by V L,n , and the normal velocity of the vapor moving away from the interface is represented by V V,n . The interface has a normal velocity of V I,n . The conservation of mass dictates that the mass flow into the interface must be equal to the mass flow out ,,,,LLnInVVnInVVVV , (2-1) where L and V are the density of the liquid and vapor phases respectively.
10 Similarly, the conservation of momentum dictates that the momentum flowing into the interface must be equal to the momentum flowing out. Since momentum is a vector quantity, the conservation of momentum in the normal direction is 22,,,,,,121122vnLnVLVLLLnInVVnInVVPPVVVVRRnn , (2-2) where is the surface tension, R 1 , R 2 are the principal radii of curvature of the interface, P V and P L are the pressures of the vapor and liquid phases, respectively. V and L represent absolute viscosities of the vapor and the liquid phases, respectively. Liquid Vv,n VI,n V L ,n s n Vapor Figure 2-4. Schematic of a liquid-vapor interface. It is assumed that the interface is very thin with constant surface tension and that the tangential velocities of each phase at the interface are equal due to the no-slip condition. Hence, the momentum conservation in the tangential direction is written as follows: ,,,LnLsVnVsLVVVVVsnsn , . (2-3) The conservation of energy condition reads as follows: ,,,,intintLLnInLVVnInVVLTVVhVVhkknn T . (2-4)
11 In this equation h L and h V are the latent heat of vaporization of the liquid and vapor phases, k L and k V are the thermal conductivities of the liquid and vapor phases, respectively. 22.214.171.124 Rayleigh-Plesset Equation Rayleigh (1917) has presented a mathematical analysis of the growth and collapse of a spherical cavity. The equation in its most general form is known as the Rayleigh-Plesset equation (Brennen 1995). A spherical bubble of radius R(t) in an infinite domain of liquid is shown in Figure 2-5. The pressure in the liquid P (t) is assumed to be known and the temperature of the liquid is assumed to be constant. The pressure P B (t) and temperature T B (t) inside the bubble are homogeneous and uniform. With spherical symmetry, the radial flow is irrotational and the velocity potential is rtRU2)(. , where d t tdR)( U. R(t) u ( r,t) Bubble Liquid P (t), T (t) Figure 2-5. Schematic of a bubble in an infinite domain of liquid. 22()()()(,) 2 R tdRtFturtrrdtr (2-5) The spherically symmetric form of the Navier-Stokes equation for a Newtonian fluid is
12 ,211222rururrrruuturPLL (2-6) where L , L are the density and the kinematic viscosity of the liquid, respectively. Substituting Eq. (2-5) into Eq. (2-6), the following is obtained 522211rFdtdFrrPL . (2-7) Integrating the above equation from r to one gets 42211rFdtdFrPPL . (2-8) The net force on the interface in the normal direction is equal to R2PBrr , where is the surface tension and rr is the normal stress in the radial direction ruPLrr 2 . (2-9) In the absence of mass transport the net force on the interface must be zero. Substituting the above equations into Eq. (2-8), one ends up with the Rayleigh-Plesset equation (Brennen 1995) 222()()3422BLLLPtPtdRdRdRRdtdtRdtR . (2-10) In the above equation, R is the radius of the spherical cavity, L is the liquid density, L is liquid kinematic viscosity, P B (t) is the pressure inside the cavity and P (t) is the pressure of the free stream. The contribution of the non-condensable gas can be incorporated into the Rayleigh-Plesset equation. In the absence of thermal effects, non-spherically symmetric internal
13 motion tends to mix the contents of the bubble. This can improve the validity of the polytropic assumption (Brennen 1995) koGoGRRPP3 , (2-11) where R o , P Go are the initial radius and pressure of the non-condensable gas, respectively. The polytropic constant k is equal to 1.0 for constant bubble temperature. Then the bubble pressure can be written as the summation of the vapor pressure and partial pressure of the non-condensable gas 3()()koBVGoRPtPTPR , (2-12) where P V (T ) is the vapor pressure at the surrounding liquid temperature. Then the Rayleigh-Plesset equation in the absence of thermal effects including the contribution of non-condensable gas becomes 3222()()342kVGooLLLpTptPRdRdRdRR 2L R dtdtRdtR . (2-13) Thermal effects can also be incorporated into the Rayleigh-Plesset equation by assuming that non-condensable gas has an ideal gas behavior. Since the mass transfer acros the interface is neglected in this model, the vapor inside the bubble can have a different vapor pressure. Initially the non-condensable gas is at radius R o and temperature T and it expands to radius R and has a temperature of T B . Then the pressure of the non-condensable gas can be written as 3RRTTPPoBGoG . (2-14)
14 2.1.5 Cavitation Inception and Cavitation Parameter A commonly used criterion for cavitation inception is based on a static approach and states that the cavities occur when the hydrodynamic pressure drops below the vapor pressure of the liquid at the free stream temperature. This can be formulated as VPP . (2-15) Joseph (1995) has proposed an improved criterion, which can account for anisotropic flow structure. The criterion is formulated in terms of the principal stresses occurring in a moving fluid rather than the pressure in a static fluid used commonly. The formulation is based on the idea that the liquid will rupture in the direction of maximum tension. It has been stated that this can lead to cavity formation under different conditions. For a Newtonian fluid the criterion can be written as follows: 2 ~ VPPxu , (2-16) where vP ~ is an empirical criterion that replaces P V to take into account the impurities for practical applications. Since the flow dynamics is taken into consideration, cavity formation can occur under different conditions than Eq. (2-15). No study has been found utilizing this criterion. A dimensionless parameter is desirable to describe the flow conditions relative to those for cavitating flows and to obtain a unique value for each set of dynamically similar cavitation conditions. In cavitation terminology, denotes the parameter of similitude ,212 VPPLV (2-17)
15 where P and V are the pressure and velocity at a reference point respectively. Unfortunately, is not strictly a parameter of similitude for cavitating flows. It is a necessary but not a sufficient condition. It is known that nuclei present in free stream flow have an impact on the two-phase flow structure (Knapp et al. 1970). 2.1.6 Kinetic Theory of Mass Transfer In order to sustain a steady cavity shape, the mass transfer into the cavity due to evaporation must be in equilibrium with the mass transfer out of the cavity due to condensation. It is also likely that mass transfer out of the cavity can occur when large segments of small cavities detach from the main cavity. Figure 2-6. Schematic of the mass transfer through a cavity (adapted from Franc et al. 1995). This process takes place in the rear part of the cavity (Franc et al. 1995). The schematic of the mass transfer through a steady cavity is given in Figure 2-6. The physics related to the interface mass transfer can be explained in the context of kinetic theory. The mass transfer through the interface is basically the arrival and departure of molecules from the interface. It is physically not correct to assume a thermostatic pressure and temperature
16 on each side of the interface because conditions close to the interface are not at static thermal equilibrium (Collier 1994). This implies that there is a temperature jump across the interface. Using kinetic theory, one can show that the flux of molecules through an imaginary plane is given by the following equation (Collier 1994) 2/12/12TPRMj , (2-18) where M is the molecular weight, R is the universal gas constant, P is the pressure and T is the absolute temperature. The net flux of molecules through an interface can be expressed as the difference of the flux of molecules due to evaporation (j ) and condensation (j + ) jjj . (2-19) Using Eq. (2-18) the net flux is 2/12/12/12LLVVTPTPRMj . (2-20) In the above equation, the subscripts L and V represents the liquid and vapor phases respectively. Assuming small differences in pressure and temperature, the above equation simplifies to the following TTPPPRTMj222/1 . (2-21) Another simplification can be made to Eq. (2-21) assuming that there is no temperature jump across the interface PRTMj2/12 . (2-22)
17 Eq. (2-22) implies that condensation occurs when the vapor pressure exceeds the liquid pressure by a certain amount. It is also known that the vapor pressure over a drop P g is higher than the vapor pressure of a planar interface P and is given by the following expression (Collier 1994) LLVvrPvPP*21 , (2-23) where r * is the drop radius, is the surface tension, v f and v g are the specific volume of the liquid and vapor phases respectively. The above discussion helps us to understand the physics of cavitation criteria. As it is described previously, cavitation occurs when the pressure of the liquid drops below the vapor pressure, which is an evaporation process. However, in order to maintain a steady cavity shape, it is also reasonable to expect a condensation process in the cavity. With the help of the kinetic theory, one can say that the condensation may occur at a pressure above the vapor pressure at which evaporation takes place. The amount of the pressure difference can be calculated from Eq. (2-22). Since condensation typically occurs at the rear portion of the cavity, the interface at this region is far from planar and this again suggests that a hysteresis may exists for the cavitation criteria. 2.2 Experimental Studies An extensive amount of literature exists on experimental studies of cavitation. Among these studies only the ones that are most recent and give insight into the flow structure of cavitation are covered in this section. Table 2-2 gives an overview of experimental studies. In the following, selected studies are summarized.
18 Arakeri and Acosta (1973) have studied the cavitation inception and development on two axisymmetric bodies using the Schlieren flow visualization method. Cavitation inception has been found to occur in the laminar boundary layer separation region. There is a close correlation between the point of laminar separation and cavity detachment. Arakeri (1975) has also shown that the distance between these two points are strongly dependent on the Reynolds number. Franc and Michel (1985) have investigated this issue further by studying the cavitating flow over circular and elliptic cylinders and NACA foils. The visualization of the boundary layer with dye injection has shown that there is a strong interaction between the attached cavitation and the boundary layer. For flows with low gas nuclei content, it has been stated that a large developed cavity always detaches downstream of a laminar separation. Unexpected results have been observed in the absence of laminar boundary layer separation, such as the disappearance of developed cavities as results of the decrease in the ambient pressure and also the disappearance of large cavities after a transition to turbulence. Kubota et al. (1989) have investigated the unsteady structure of cloud cavitation using a conditional sampling technique. The unsteady flow velocity around the cloud cavitation was measured by a Laser Doppler Anemometry. The velocity field has been divided into large-scale and small-scale components. The large scale unsteady cloud cavitation motion has been associated with the large-scale velocity. The small-scale velocity has been found to be active on the boundaries of the cloud cavitation. Furthermore the cloud cavitation moves with a velocity smaller than the mean free stream velocity and it has a concentrated vorticity region at its center.
19 Kuhn de Chizelle et al. (1995) have studied the traveling bubble cavitation and its interactions with the flow. Cavitating flows over Schiebe headforms have been investigated using high-speed film and electrode sensors to measure noise signals. It has been found that the micro-fluid mechanics is important and the spherical Rayleigh-Plesset analysis may be insufficient to describe the process. Micro fluid mechanics phenomena involves thin films separating the bubbles from the surface, formation of attached streaks in the wake of the bubbles and the process of bubble fission and roll-up during the collapse of the bubble. The measured impulses during bubble collapse have shown that the noise levels are lower than the levels predicted by the Rayleigh-Plesset analysis. It has been argued that the shearing and the fission of the bubbles are responsible for this difference. Li and Ceccio (1996) have studied along the same lines and investigated the interaction of traveling bubbles with the boundary layer using the Particle imaging velocimetry. It has been shown that traveling bubbles close to the surface can cause local turbulent regions by squeezing the boundary layer and can create streamwise vorticity. These local regions of turbulent flow can also sweep away a portion of the attached cavity. Kawanami et al. (1997) have conducted a series of experiments to control the cloud cavitation. It has been demonstrated that re-entrant jet downstream of the cavity is the basic mechanism that triggers the shedding of the cloud cavitation. This process can be controlled by adding a small obstacle to block the re-entrant jet. Cavitation noise and drag can be decreased effectively by this modification on the hydrofoil surface. The velocity of the re-entrant jet is estimated as to be the same order as the magnitude of free stream flow.
20 Ceccio et al. (1997) have performed experiments to understand the effect of salt and fresh water on cavitation inception, bubble size and acoustic emission. For identical flow conditions, significant differences have been found in fresh and salt water. The bubble size and event rate of cavitation are reduced, on the other hand acoustic emission increases under salt-water conditions. The research has shown that the chemical differences and the free stream nuclei distribution has an impact on cavitation process and the authors hypothesize that the salt solution reduces the number and size of free stream nuclei. Stutz and Reboud (1997a) have studied the two-phase flow structure of sheet cavitation. The void fraction, flow velocity and bubble size inside the cavities have been measured. Interaction of cavities with the outer flow has also been investigated. The Reynolds number is in the range of 4x10 5 -2x10 6 and the cavitation number is in the range of 0.60 to 0.75. The measurements have indicated that pressure inside the cavity is close to vapor pressure of the liquid at the corresponding temperature and the pressure downstream of the cavity. The streamlines around the cavity indicates that the outer flow passes through the interface and the mass flow rate within the cavity reaches 4% of the total mass flow rate. Even though there is no reversed flow under non-cavitating conditions, it is present in the downstream part of the sheet cavitation. Reisman et al. (1998) have investigated the acoustics of cloud cavitation using piezo-electric pressure transducers and high-speed video camera. Several types of propagating structures, so called bubbly shock waves, have been observed. The largest impulsive and radiated noise has been attributed to the coherent collapse of a well-defined separate cloud in regions of high pressure.
21 Gopalan and Katz (2000) have looked into the flow structure in the closure region of sheet cavitation. It has been demonstrated that the collapse of vapor bubbles is the primary source of vorticity production. It has been found that small changes in the size of cavity (due to changes in cavitation parameter) result in substantial changes in turbulence level and momentum thickness of the downstream boundary layer. The results obtained are also used to check the mass and momentum balance downstream of the cavity. Callenaere et al. (2001) have studied the cavitation instability induced by the development of a re-entrant jet. Two main classes of instabilities have been proposed: Intrinsic instabilities and system instabilities. The intrinsic instability is the one that originates within the cavity itself, whereas the system instability originates due the interaction of the cavity with other parts of the hydraulic system. It has been found that if the adverse pressure gradient at the closure region of the cavity is strong, then a re-entrant jet develops. The cavity behavior also depends on the relative thicknesses of the cavity and the re-entrant jet. The interaction of the re-entrant jet with thick enough cavities leads to cloud cavitation. If the relative thickness of the cavity is comparable to re-entrant jet thickness then the cavity appears as two-phase mixture.
22 Table 2-2. An overview of the past experimental studies on cavitating flows. Reference Experimental method Problem Studied Main Findings Arakeri and Acosta (1973) Schlieren flow visualization method Cavitating flow over two axisymmetric bodies Close correlation between the point of laminar separation and cavity detachment. Arakeri (1975) Schlieren flow visualization method Cavitating flow over two axisymmetric bodies Distance between point of laminar separation and cavity detachment is strongly dependent on Reynolds number. Franc and Michel (1985) Flow visualization by dye injection Cavitating flow over circular and elliptic cylinders and NACA foils Strong interaction between attached cavitation and boundary layer. Cavity detaches behind a laminar separation. Kubota (1989) LDV and conditional sampling technique Cavitating flow over hydrofoils Large scale unsteady cloud cavitation is associated with large scale velocity. Concentrated vorticity region at the center of cloud cavitation. Kuhn de Chizelle (1995) High-speed films and electrode sensors for noise signals Cavitating flow over Schiebe headforms Rayleigh-Plesset analysis may be insufficient to describe the micro-fluid mechanics. Observations of bubble fision and roll up during collapse. Li and Ceccio (1996) Particle imaging velocimetry Cavitating flow over 2-D hydrofoils Traveling bubbles close to the surface can cause local turbulent regions. Kawanami et al. (1997) High-speed video and photo, pressure measurements by pressure pick ups and hydrophone Cloud cavitation on a hydrofoil section A re-entrant jet from downstream to upstream triggers the collapse of sheet cavitation. It can be controlled by adding a small obstacle on the hydrofoil. Ceccio et al. (1997) Electrodes, piezoelectric hydrophone Cavitating flow over axisymmetric head forms For identical flow conditions, bubble cavitation has different characteristics in fresh and salt water. Stutz and Reboud (1997a) LDV and double sensor optical probe technique Two-phase flow structure of sheet cavitation on venturi type test section Measured the distribution of mass and momentum balance inside the cavity and also the mass flow rate within the cavity. Reisman et al. (1998) Electric pressure transducers and high-speed video camera Cloud cavitation over hydrofoils Observations of bubbly shocks during collapse of the cavitation cloud. Gopalan and Katz (2000) Particle image velocimetry and high speed photography Flow structure in the closure region of attached cavitation in a nozzle. Collapse of cavities in the closure region results in vorticity production. Highest turbulence levels are located downstream of the closure region. Callenaere et al. (2001) LDV measurements Cavitation behind a diverging step in a venturi test section. Two types of cavitation instability: system and intrinsic. Strong adverse pressure gradients are required for re-entrant jet development.
23 2.3 Numerical Studies Computational modeling of cavitation has been pursued for years. Early studies primarily utilize the potential flow theory; they are still widely used in many engineering applications. Studies dealing with cavitation modeling through the computation of the Navier-Stokes (N-S) equations have emerged in the last decade. An overview is given in Table 2-3. These studies can be put into two categories, namely interface tracking methods and homogeneous equilibrium flow models. In the first category the cavity region is assumed to have a constant pressure equal to the vapor pressure of the corresponding liquid and the computations are done only for the liquid phase. Constant pressure assumption is physically sensible and has been verified experimentally (Shen and Dimotakis 1989, Stutz and Reboud 1997a). Computationally, the liquid-vapor interface can be tracked based on this assumption, along with a wake model to handle the shape of the cavity. Grid is often regenerated iteratively to conform to the cavity shape. These models are capable of simulating sheet cavitation but may not be adequate for cases where bubble growth and detachment exists. In addition, so far, they are limited to 2-D planar or axisymmetric flows because of the difficulties involved in tracking 3-D interfaces. Examples can be found in Chen and Heister (1994), and Deshpande et al. (1997). The second category can be termed the homogeneous equilibrium flow models in which the single-fluid modeling approach is employed for both phases. Differences between the various models in this category mostly come from the relation that defines the variable density field. Delannoy and Kueny (1990) utilize an arbitrary barotropic equation of state to compute the density field. Likewise, Chen and Heister (1996) derive a time and pressure dependent differential equation for density. Ventikos and Tzabiras
24 (2000) introduce the water-vapor state laws to model the cavitation dynamics, and consider the whole domain, including both vapor and liquid phases, as a compressible fluid. A pressure-based algorithm is adopted, where the density correction in vapor phase is linked to the mass residual there. The results obtained do not seem to recover the cavity vapor pressure in accordance with the cavitation number. Edwards et al. (2000) use the Sanchez-Lacombe equation of state and solve the temperature transport equation in addition to the Navier-Stokes equations. Considering the isothermal character of cavitating flows in many applications, utilizing a temperature or enthalpy equation and assuming the whole flow field as compressible may not be the most effective approach. Kubota et al. (1992) couple the Rayleigh-Plesset equation to the flow solver and compute the void fraction based on the bubble radius. Then the density is calculated using the void fraction formalism. Due to the time dependent nature of the Rayleigh-Plesset equation, the model is restricted to unsteady cloud cavitation. The authors have also reported that this method is prone to instability because of high pressure-density dependence, and could not reach the convergence levels of noncavitating flow simulations. To account for the cavitation dynamics in a more flexible manner, recently, a transport equation model is developed. In this approach volume or mass fraction of liquid (and vapor) phase is convected. Singhal et al. (1997), Merkle et al. (1998), and Kunz et al. (2000) have employed similar models based on this concept with differences in the source terms. One apparent advantage of this model comes from the convective character of the equation, which allows modeling of the impact of inertial forces on cavities like elongation, detachment and drift of bubbles. Singhal et al. (1997) utilize a pressure-based algorithm but offer no detailed information related to computational convergence and
25 stability. Merkle et al. (1998) and Kunz et al. (2000) have employed the artificial compressibility method with special attention given to the preconditioning formulation. Ahuja et al. (2000) have developed an algorithm to account for the compressibility effects in the context of artificial compressibility methods with adaptive unstructured grids. Venkateswaran et al. (2001) have compared the above-mentioned studies, and concluded that all three preconditioning formulations are essentially the same with minor differences.
26 Table 2-3. An overview of the past numerical studies on cavitating flows. Reference Cavitation model Numerical algorithm Applications/Main Findings Kubota et al. (1992) Rayleigh-Plesset equation coupled to the Poisson equation. Cavity region is modeled as compressible fluid with variable density Marker and Cell 3-D N-S equations No turbulence model. Cloud cavitation on hydrofoils. Numerical instability for high-density ratio. Re=3x10 5 . Chen and Heister (1994) Interface tracking based on P=Pvap Grid conforms to the cavity shape. Marker and Cell 2-D N-S equations No turbulence model Pressure distribution on axisymmetric geometries. Re=1.36x10 5 . Chen and Heister (1996) Time and pressure dependent pseudo-density equation. Marker and Cell 2-D N-S equations No turbulence model Pressure distribution on axisymmetric geometries. Re=1.36x10 5 . Deshpande et al. (1997) Interface tracking based on P=Pvap with mass transfer. Grid conforms to the cavity shape. Artificial Compressibility 2-D N-S equations No turbulence model Sheet cavitation for cryogenic fluids. Studied the thermal boundary layer over cavity. Singhal et al. (1997) Vapor mass fraction equation with pressure dependent source terms. Pressure-based 2-D N-S equations kturbulence model Pressure distribution and discharge coefficient for orifices and hydrofoils. Re=2.x10 6 . Merkle et al. (1998) Vapor mass fraction equation with pressure dependent source terms. Artificial Compressibility 2-D N-S equations Two equation turbulence model Pressure distribution on hydrofoils. Kunz et al. (2000) Volume fraction equation with pressure dependent source terms Nonconservative continuity equation. Preconditioning strategy. Artificial Compressibility 3-D N-S equations kturbulence model Pressure distribution on axisymmetric geometries. Re=1.36x10 5 . Ahuja et al. (2000) Vapor mass fraction equation with pressure dependent source terms. Preconditioning strategy. Adaptive unstructured meshes Artificial Compressibility 3-D N-S equations kturbulence model Simulations of cavitating flow over hydrofoils (Re=2x10 6 ) and axisymmetric geometries. (Re=1.36x10 5 ). Edwards et al. (2000) Temperature distribution is computed to determine density variation based. Sanchez-Lacombe equation of state. Artificial Compressibility 3-D N-S equations Spalart-Allmaras one-equation model Pressure distribution on axisymmetric geometries. Reported poor convergence and pressure overshoots Re=1.36x10 5 . Ventikos and Tzabiras (2000) Temperature distribution is computed to determine density variation based on steam-water tables. Pressure-based 2-D N-S equations No turbulence model Pressure distribution over airfoils. Re=2000 in computations while Re=2.5x10 6 in experiments. Venkateswaran et al. (2001) Discussed the preconditioning strategies utilized in Kunz et al. and Ahuja et al. Artificial Compressibility 3-D N-S equations kturbulence model Pressure distribution on axisymmetric geometries. Re=1.36x10 5 .
CHAPTER 3 THEORETICAL FORMULATION In this chapter the governing equations and the computational methodology are presented. First, the Navier-Stokes equations governing the flow of incompressible single-phase fluids are presented. Second, the Favre-averaged Navier-Stokes equations and the turbulence modeling concept are provided. Following these, cavitation modeling is discussed, and a transport equation-based interfacial dynamics cavitation model is developed. Finally, the chapter is concluded with the presentation of numerical methodology and the pressure-based algorithms for steady and time-dependent computations. The governing equations are presented in Cartesian coordinates for the ease of discussion. Issues related to the implementation of the models are also discussed. 3.1 Governing Equations The Navier-Stokes equations governing the flow of a Newtonian fluid in the absence of body forces are presented below: 0)( ut (3-1) ijPuuut )()( , (3-2) where viscous stress tensor is given by the following equation kkijijjiijxuxuxu32 . (3-3) 27
28 The above equations derive from the physical conservation laws of mass and momentum. For a complete derivation of these partial differential equations, one can refer to any fluid mechanics textbook (e.g., Johnson 1998). Understanding of fluid turbulence is not complete and the prediction of turbulence continues to advance. The above equations are general and can, in principal, be applied to do a direct numerical simulation (DNS) of turbulent flows. Such a task requires that all relevant scales of the turbulent flow be resolved; hence, great amount of computer resources are needed. To date, DNS of only simple flows at modest Reynolds numbers have been performed and DNS of complex high Reynolds number flows are beyond the capabilities of current computing resources (Pope 2000). On the other hand, most of the time, the interest is on the mean flow variables. Once the attention is focused on the mean quantities, then it is more convenient to work with the mean flow equations. In Reynolds averaging, the instantaneous flow variables are decomposed into a mean ( ) and a fluctuating component ( ). For general turbulence, the mean of a quantity is defined by ensemble averaging as follows: ),(1lim1)(txNNnnN . (3-4) The averaging is done for N realizations of the flow field. The mean of the fluctuating component is by definition equal to zero. For variable density flows, Reynolds averaging leads to additional terms involving the products of fluctuations between density and other variables. To avoid this, density weighted averaging (Favre-averaging) is introduced. Density-weighted average is defined as follows: ~ , . (3-5) ~
29 It is important to note that the mean of the double-primed fluctuation is not equal to zero but instead the mean of the double-primed fluctuation times the density ( ) is zero. Further details of the averaging process are given in Tannehill et al. (1997) and Johnson (1998). Having made these definitions, the Favre-averaged Navier-Stokes equations read the following 0)~( ut (3-6) )()~~()~(jiijuuPuuut (3-7) kkijijjikkijijjiijxuxuxuxuxuxu32~32~~ . (3-8) Note that the fluctuations of viscosity are neglected in writing the ensemble-averaged viscous stress tensor. In practice, the second term that appears in ensemble-averaged viscous stress tensor is also neglected on the basis of order of magnitude analysis (Tannehill et al. 1997). It should be mentioned that during the averaging process some of the information is lost and because of the nonlinear structure of the N-S equations additional unknowns namely Reynolds stress terms ( jiuu ) arise. As a result, the set of equations are not closed due to the additional six unknowns in this symmetric tensor. In order to close the system of equations, a modeling concept needs to be provided to relate the Reynolds stresses to the mean velocity field. The Boussinesqâ€™s eddy-viscosity hypothesis is one of the oldest proposals for the turbulence closure (Schlichting 1979). It can be put in the following form kxuxuxuuukktijijjitjiRij~~32~~ . (3-9)
30 Utilizing the Boussinesqâ€™s eddy-viscosity hypothesis the Favre-averaged Navier-Stokes equations read the following 0)~( ut (3-10) )()~~()~(RijijPuuut (3-11) kkijijjitRijijxuxuxu~32 ~ ~)( . (3-12) Note that the structure of the above Favre-averaged Navier-Stokes equations are the same as the Reynolds-averaged Navier-Stokes equations and as well as the original Navier-Stokes equations, and hence, makes it convenient for computational purposes. In the present Navier-Stokes solver, the governing equations written in body-fitted curvilinear coordinates are solved. The transformation of governing equations into body-fitted curvilinear coordinates is well established and the related information is documented in details in Shyy (1994) and Thakur et al. (1997). 3.2 Turbulence Modeling After obtaining the Reynolds or Favre averaged Navier-Stokes equations the task is to estimate the eddy-viscosity for which a turbulence model has to be incorporated into the course of computation. A review for the assessment of turbulence modeling is given by Spalart (2000). Commonly known models can be classified as follows: Zero and one equation models Two equation model; kmodel, kmodel Reynolds stress model Algebraic stress model Large eddy simulation
31 3.2.1 Two Equation (k-) Turbulence Model Among the mentioned models, kmodel (Jones and Launder 1972) has become the popular one with well-known deficiencies. It is computationally tractable and robust. A wide variety of flow structures have been computed with this model (Wilcox 1993). In this model, two partial differential equations, one for the turbulent kinetic energy (k) and a second for the rate of dissipation of the turbulent kinetic energy () are solved, where k and are defined as follows (Wilcox 1993): iiuuk21 jiijxu . (3-13) The transport equation for these properties is presented below: kuuktkktRij)(~)()~( (3-14) )()~(221tkCPkCut (3-15) uPRij ~ )( 2kCt . (3-16) In the above equations C 1 , C 2 , k , are empirical constants. For certain types of flows, where flow equilibrium does not exist such as flows with recirculation, rotation and large streamline curvatures, these constants have to be modified to account for the nonequilibrium effects. The coefficients C 1 and C 2 regulate the production and dissipation in the equation and the magnitude of eddy viscosity, respectively. There are other variants of the original kturbulence model introduced by modifying the empirical constants. Only one of them along with the original model is considered in this study. The adopted models are designated as original kand nonequilibrium kfor the ease of
32 study. The empirical constants used in these models are tabulated in Table 3-1. Further information about these models and validation for different test cases are given in Shyy et al. (1997). Table 3-1. The empirical constants used in the variants of kturbulence model. P is the first term on the right hand side of Eq. (3-15). MODEL C C 1 C 2 k Original k0.09 1.44 1.92 1.0 1.3 Nonequilibrium k0.09 1.15+0.25(P/) 1.45+0.45(P/) 0.8927 1.15 3.2.2 Wall Functions The implementation of kturbulence model requires a special attention in the vicinity of a solid wall. Physically, the flow near the no-slip boundaries is different than the free stream flow. Viscous effects dominate in the near-wall region. Because of this difference in the flow structure, a treatment is needed along with the kturbulence model. A common approach to resolve the near-wall turbulence is the so-called wall functions (Launder and Spalding 1974). The nondimensional normal distance from the wall, y + (local Reynolds number), and the nondimensional velocity, u + , can be used to correlate the profiles in the near-wall region of turbulent flows *uuut wallu* *uyy , (3-17) where u * is the friction velocity and is the kinematic viscosity. A detailed discussion on the derivation of these terms is given in Schetz (1984). The above definition of y + is referred to as â€œoriginal definitionâ€. The wall function treatment for fully developed turbulent flow in the vicinity of a solid wall, assuming that flow has two-layer structure (viscous sublayer followed by the log layer), is as follows (White 1974):
33 63.11)log(163.11yEyyyu , (3-18) where E and are constants. The above definitions for the nondimensional quantities cause a singularity when there is separation in the flow. To eliminate the computational singularity, one can use the concept of local equilibrium between production and dissipation of turbulent kinetic energy. The following can be written (Shyy et al. 1997): kCwall . (3-19) This formula relates the wall shear stress to the turbulent kinetic energy. Since k is not zero at separation, the above form does not cause computational singularities. Using the flow equilibrium assumption, one can define the nondimensional normal distance, referred to as â€œequilibrium definitionâ€, as follows: PykCy4/1 , (3-20) where y P is the normal distance of the first grid point away from the wall. The wall shear stress components in the u-, v-, w-momentum equations are as follows: txPxwalluuyy tyPywalluuyy tzPzwalluuyy . (3-21) A detailed discussion of the wall functions and its implementation in body fitted coordinates is provided in Shyy et al. (1997). 3.2.3 Implications of the Near-wall Turbulence Treatment In this section, the outcomes of different y + definitions plus the equilibrium of the flow structure and the impact of grid distribution are demonstrated. Turbulent flow calculations are performed for a flat plate. Two different grids are tested for angles of
34 attack of 0 and 5, respectively. The Reynolds number based on the length of the plate is 10 6 . Figure 3-1 shows the grid distributions for the first 11 grid points in j direction. In the dimensions and some features of the grids are summarized. Figure 3-2 shows the wall shear stress profiles obtained from computations on two different grids. The computational results are compared with the following correlation (Schlichting 1979) Table 3-2 Table 3-2. Description of grid properties. 512)(Re0148.0xwallU . (3-22) Grid dimensions (i,j) Number of grid points on the wall Grid-A 191x81 51 Grid-B 191x81 71 The wall shear stress profile of the case with Grid B is closer to the empirical result, because it has a smoother grid distribution, more grid points, and has the first grid point above the flat plate in the log layer. Hence, even though the wall function is designed to match the velocity profile in a range of the y + values, the accuracy of the computation is noticeably affected by the distributions of the grid and y + values in the near-wall region. In Figure 3-3 and Figure 3-4, one can see that different definitions of y + give almost the same result. The reason for this is the equilibrium between the production and dissipation of turbulent kinetic energy on the flat plate. However, y + values are not the same when there is an angle of attack. Since the flow is not in equilibrium, y + values at the leading edge of the flat plate, as indicated in Figure 3-5, are far from being close. The difference diminishes at the trailing edge, which is an indication of flow equilibrium. Even if there is a significant difference in y + profiles; the effect on wall shear stress profile is not large as shown in . Thus it can be concluded that y + has a secondary effect on wall shear stress. This is also evident from the wall shear stress terms of momentum equation, Eq. (3-21). The wall shear stress is proportional to the ratio Figure 3-6
35 of y + to u + . Since u + is a function of y + the overall effect on wall shear stress is not dominant. Grid A (191x81) Grid B (191x81) 1 2 3 4 5 6 7 8 9 10 11 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 y coordinate j index Figure 3-1. Grid distribution along j direction. Figure 3-2. The effect of grid distribution on wall shear stress.
36 Figure 3-3. y + profiles along the flat plate. =0, Grid B is used. The kturbulence model is popular to predict fully turbulent flows because it is computationally economic. Its capabilities and deficiencies are also well documented. The goal of in this section is to assess the predictive capability of the kmodel for attached and mildly separated flows. In the kturbulence model, the wall function treatment on different grids with varying angles of attack, between 0 and 5 is investigated. Issues related to the implementation of the equilibrium approximation and the grid distribution are important because they noticeably affect the computational accuracy.
37 Figure 3-4. Effect of different y + definitions on wall shear stress. =0, Grid B is used. Figure 3-5. y + profiles along the flat plate. =5, Grid B is used.
38 Figure 3-6. Effect of different y + definitions on wall shear stress. =5, Grid B is used. 3.3 Cavitation Modeling A common approach in cavitation modeling is to use the homogeneous flow theory. In this theory, the mixture density concept is introduced and a single set of mass and momentum equations are solved. Different ideas have been proposed to generate the variable density field; a review of such studies is discussed in Section 2.3. Some of the existing studies solve the energy equation and determine the density through suitable equations of state (Edwards et al. 2000). Since most cavitating flows are isothermal, arbitrary barotropic equations have been proposed to supplement the energy consideration (Delannoy and Kueny 1990). Another popular approach is the transport equation-based model (TEM) (Singhal et al. 1997, Merkle et al 1998, Kunz et al. 2000, Senocak and Shyy 2001). In TEM, a transport equation for either mass or volume fraction, with appropriate source terms to regulate the mass transfer between phases, is solved. Different source terms have been proposed by different researchers, which will be
39 discussed in the next section. A recent experimental finding (Gopalan and Katz 2000) helps assess the adequacy of the above-mentioned physical models. It shows that vorticity production is an important aspect of cavitating flows, especially in the closure region (Gopalan and Katz 2000). Specifically, this vorticity production is a consequence of the baroclinic generation term of vorticity equation P1 . (3-23) Clearly, if an arbitrary barotropic equation )(Pf is used, then the gradients of density and pressure are always parallel; hence the baroclinic torque is zero. This suggests that physical models that utilize a barotropic equation will fail to capture an experimentally observed characteristic of cavitating flows. Likewise, solving an energy equation will also experience the same situation if the flow is essentially isothermal. However, in TEM approach, the density is a function of the transport process. Consequently, gradients of density and pressure are not necessarily parallel, suggesting that TEM can accommodate the baroclinic vorticity generation. 3.3.1 Transport Equation-Based Empirical Cavitation Models Different forms of the transport equation-based cavitation models have been proposed in literature. The main differences are due to different source terms that are needed for cavity generation and destruction. All the models have introduced empirical factors to regulate the mass transfer. These empirical factors have been tuned based on numerical experimentation. So far, satisfactory results have been produced with these models, provided that the computational algorithm addresses the numerical issues specific to cavitating flows. However, the derivation of these models has not been clarified in the literature and no satisfactory physical arguments have been made regarding the source of
40 the empirical factors. In what follows, these models are presented and discussed. In Section 4.2, the predictive capability of these models is assessed. 126.96.36.199 Model-1 (Singhal et al. 1997, Merkle et al. 1998) Several researchers have adopted this model (Singhal et al. 1997, Merkle et al. 1998, Ahuja et al. 2001). Both volume fraction and mass fraction forms of it have been adopted. Evaporation and condensation terms are both functions of pressure. The liquid volume fraction form is considered in this study. tUPPMAXCtUPPMINCutLLVLprodLVLVLLdestLL)50.0(10,)50.0(0,)(22 (3-24) The empirical factors adopted in this study have the following values, C dest =1.0, C prod =8.0x10 1 . These parameters are determined through numerical experimentation and have been non-dimensionalized with free stream values to have the correct dimensional form as the convective terms. The same parameters can be used for different geometries and flow conditions provided that they are non-dimensionalized with the free stream values (Senocak and Shyy 2001). The discussion holds for the other cavitation models also. 188.8.131.52 Model-2 (Kunz et al. 2000) The liquid volume fraction is chosen as the dependent variable in the transport equation. The evaporation term is a function of pressure whereas the condensation is a function of the volume fraction. tCtUPPMINCutLLLprodLLLVLVdestLL1)50.0(0,)(22 . (3-25) The empirical factors adopted in this study have the following values, C dest =9.0x10 5 , C prod =3.0x10 4 .
41 184.108.40.206 Model-3 (Singhal et al. 2001) The vapor mass fraction is the dependent variable in the transport equation. Both evaporation and condensation terms are functions of pressure. The model equations that are presented here are slightly different from the ones in the original paper (Singhal 2001), which considers non-condensable gas. No such limitation is imposed here. )()(mmuftfVmVm (3-26) VLVVLprodVLVVLdestPPifPPUCmPPifPPUCm:32:322121 . (3-27) The empirical factors adopted in this study have the following values, C dest /=1.225x10 3 , C prod /=3.675x10 3 . 3.3.2 Derivation of a Transport Equation-Based Interfacial Dynamics Cavitation Model The derivation starts by considering a liquid-vapor interface. The mass and normal momentum conservation at such an interface can be written as follows (Carey 1992): )()(,,,,nInVVnInLLVVVV (3-28) 2,,2,,,,21)()(2211nInVVnInLLnLLnVVLVVVVVnVnVRRPP (3-29) Figure 3-7 illustrates a typical representation of a liquid-vapor interface based on the homogeneous flow theory. The mixture density is defined as follows based on the liquid volume fraction (Wallis 1969) )1(LVLLm . (3-30)
42 U L V m Figure 3-7. Representation of a vaporous cavity in homogeneous flow theory As shown in Figure 3-7, a hypothetical interface is assumed to lie in the liquid-vapor mixture region. One can further assume that the phase change takes place between the vapor and the mixture phase, because in the sharp interface case, Eq. () can represent either the vapor or liquid density. In this formulation, the liquid phase is represented by the mixture phase and hence, the mass and normal momentum conservation equations reduce to the following forms after neglecting the viscous and surface tension effects 3-30 3-30 )()(,,,,nInmmnInVVVVVV (3-31) 3-31 2,,2,,)()(nInVVnInmmLVVVVVPP . (3-32) From the mass conservation condition, Eq. (), the following relation can be deduced mnInVVnInmVVVV )()(,,,, . (3-33) The momentum conservation condition, Eq. (3-32), can be rearranged to the following form by incorporating the mass conservation condition, Eq. (3-31), and Eq. (3-33). 1)(2,,mVnInVVLVVVPP . (3-34) At this point the definition of mixture density, given in Eq. (), is incorporated into the above equation that leads to the following forms
43 1)1()(2,,LVLLVnInVVLVVVPP (3-35) 2,,)()1()()()(nInVVLVLVLLLVLLVVVPPPP . (3-36) The final form reads the following after further arrangement of the terms )()()1)(()()()(2,,2,,VLnInVLVLVLnInVVLVLLLVVPPVVPP . (3-37) The above equation defines the value of liquid volume fraction. In the context of TEM, it is necessary to couple the above interfacial condition as a source term into the transport equation of liquid volume fraction. Many researchers in turbulence modeling have used the dimensional analysis and physical reasoning as closure approximations (Wilcox 1993). The same approach is adopted here and Eq. (3-37) is simply normalized with a characteristic time scale defined as ULtch so that a dimensionally correct form that represents the rate of generation of L appears tVVPPtVVPPtSVLnInVLVLVLnInVVLVLLL)()()1)(()()()(2,,2,, . (3-38) The above source term coupled to the transport equation of L is shown below tVVPPtVVPPutVLnInVLVLVLnInVVLVLLLL)()()1)(()()()()(2,,2,, . (3-39) No empirical constants appear in the above equation. The first term on the right hand side, compared to the second term, is scaled naturally by a factor of the nominal density ratio ( L / V ). To utilize the above equation, the vapor phase velocity normal to the interface (V V,n ) and the velocity of the interface (V I,n ) are needed.
44 The derivation of the model is based on an existing interface; hence, conditional statements are required on the pressure terms in order to couple the model to the flow computation. The cavitation inception condition suggests that cavitation incepts once the hydrodynamic pressure drops below the thermodynamic vapor-pressure value of the corresponding liquid. As seen from Eq. (3-38), in the pure liquid phase ( L =1) the second term will be zero and only the first term can respond to a pressure drop below the thermodynamic vapor-pressure value. Hence, the inception condition is imposed as minimum (MIN) function on pressure difference term of the first term of Eq. (3-38). Then similarly, in the pure vapor phase ( L =0) the first term will be zero and only the second term can respond to a pressure increase above the thermodynamic vapor-pressure value. Hence, the desinent cavitation condition is imposed as a maximum (MAX) function on the pressure difference term of the second term. As a result the first term of Eq. (3-38) is responsible for conversion from liquid phase to vapor phase (evaporation), and the second term of Eq. (3-38) is responsible for conversion from vapor phase back to liquid phase (condensation). With these inclusions, the model equation to be solved along with the NavierStokes equations is the following: tVVPPMAXtVVPPMINutVLnInVLVLVLnInVVLVLLLL)()()1)(0,()()()0,()(2,,2,, . (3-40) Eq. (3-40) forms a transport equation-based cavitation model for Navier-Stokes computations of cavitating flows. Unlike the existing transport equation-based models in the literature, the current model does not require empirical constants to regulate the mass transfer process. The mass transfer is regulated based on the interfacial dynamics. Indeed, a careful examination reveals that the new model has similar terms as the Model-1, but the empirical factors are now replaced by the following terms
45 2,,22,,2))((150.0,))((150.0nInVVLLprodnInVVLLdestVVUCVVUC . (3-41) It is general to handle the time-dependency since the interface velocity (V I,n ) is taken into account in the model. For time-dependent computations, the model requires that an interface be constructed in order to compute the interface velocity (V I,n ), as well as the normal velocity of the vapor phase. However, in sheet type of attached cavitation the cavity is often modeled in a steady-state computation. Hence, the interface velocity (V I,n ) is zero for steady-state computations and the normal velocity of the vapor phase can be computed by taking the gradient of the liquid volume fraction (Brackbill et al. 1992, Francois 1998). The vapor phase normal velocity is the dot product of the velocity and the normal vector LLn V nunV , . (3-42) Since a body-fitted curvilinear coordinates are adopted in the computational algorithm, depending on the geometry, the grid may not be Cartesian. Hence, the normal vector of the interface needs to be computed based on the body-fitted curvilinear coordinates as shown below: xxxxLLLL (3-43) yyyyLLLL (3-44) zzzzLLLL . (3-45) The derivatives of with respect to body-fitted curvilinear coordinates are computed based on central differencing of the neighboring cell-centered nodes
46 3332312322211312111fffffffffJzyxzyxzyx (3-46) zyxzyxzyxzyxzyxzyxJ (3-47) yxyxfzxxzfyzzyfyxxyfxzzxfzyyzfxyyxfzxxzfyzzyf 333231232221131211 . (3-48) Since the interfacial velocities appear in the denominator of the source terms, there is a possibility of division by zero or by very small numbers. To avoid this, a simple filter is constructed for smoothing the interfacial terms. In the filter, first, a cutoff is applied on the interfacial term and then interfacial term at a node is updated by the average of the surrounding nodes. As mentioned above, the interface velocity is needed for time-dependent problems and it requires additional methods to track the movement of the interface. Development of such an interface tracking scheme for turbulent multi-phase flows in body-fitted curvilinear coordinates deserves a separate study and it is beyond the scope of this study. However, a realistic value can suffice for the current purposes and it can be derived from the mass conservation condition, given in Eq. (3-28), as follows: VLnLVLnVnIVVV 1,,, . (3-49) In this study, a simplification is adopted by introducing a single relationship between the liquid-phase normal velocity (V L,n ) and the vapor-phase normal velocity (V V,n ), which
47 helps supply the interface velocity information needed in Eq. (3-40) for time-dependent computations. With this simplification, Eq. (3-49) reads the following: nVVLVLnInVnLVfVVfV,,,,11. . (3-50) A value of -0.90 is used for f. 3.4 Numerical Methodology The baseline Navier-Stokes solver, documented in Shyy (1994), Shyy et al. (1997) and Thakur et al. (1997), employs a pressure-based algorithm and a finite volume approach to solve the fluid flow and energy equations, written in body fitted curvilinear coordinates, on collocated, multi-block grids in 2D and 3D domains. For the present study, the transport equation-based cavitation models, described earlier, are implemented into the solver and related modifications, regarding the convection schemes and the pressure-based algorithm, have been made for both steady and time-dependent computations. To help describe the underlying algorithm, the steady-state generic transport equation is adopted in vector form as qu )()( , (3-51) where is the generalized dependent variable, is the diffusion coefficient, and the second term on the right hand side represents the source term for the transported quantity . The above equation is transformed to an integral form, suitable for finite volume discretization, using the divergence theorem, VdqdSnuvs ).( , (3-52) which upon integration yields the following equation
48 jijijijijijibFFFFF,21,21,,21,21,][ , (3-53) where b i is the integrated form of the source term, and F represents the flux of at each control volume face and is composed of a convective and a diffusive part as follows: cfdiffcfconvcfcfSnuFFF ).( . (3-54) The diffusive flux is discretized using the second order central difference scheme, whereas the choice of discretization scheme for the convective flux often depends on the flow conditions and fluid physics (Shyy 1994). In this study two schemes are considered, namely the first order upwind (FOU) and the second order controlled variation scheme (CVS) (Shyy and Thakur 1994a, Shyy and Thakur 1994b). The first order upwind scheme is employed mainly to probe the sensitivity of the solution with respect to the choice of the convection operators. The CVS scheme is employed in all cases for experimental validation. Both convection schemes are illustrated considering the one-dimensional transport equation of the variable. In the first order upwind (FOU) scheme the value of the dependent variable is estimated using the upwind neighbor value. If one lets f i+1/2 be the first order flux at a control volume face, determined through first order extrapolations of two immediate neighboring cells, then the scheme can simply be written as )0,()0,(2121112121121niininiiniiuMaxuMaxf , (3-55) where u is the velocity component in the corresponding direction. Higher order spatial accuracy can be obtained by employing more grid points for extrapolation. However, it is known that second or higher order accurate schemes for convection terms produce oscillations around discontinuities. Later, in the results section,
49 it will be shown that shock-like discontinuities in density profiles do appear in cavitating flows. If a linear second order upwind scheme is used, oscillations in the vicinity of discontinuities can interfere with the mass transfer cavitation model and lead to unrealistic solution or even divergence. For this reason, the numerical methodology must introduce a convection scheme that does not generate nonphysical oscillations in the vicinity of sharp gradients and discontinuities. The second order accurate CVS scheme is based on the total variation-diminishing (TVD) concept (Harten 1984) and is suitable for the present problem. In the CVS scheme, the convective flux is estimated using the second-order TVD scheme to improve the formal order of accuracy and the local characteristic speeds are assigned according to the value of the local convective speed. The second order net flux term for the linearized implicit version of the CVS is presented as follows: niiiiiiniiiiiiniiniiiiniiniiiiiiQbrQbrrQbrQbff)()(41)()(41)()(21121)()(211212123232323122323232311212121211121212121)2(21)2(21 . (3-56) The superscript (2) stands for second order accuracy, (n+1) represents the value at the current iteration, and subscript (i) indicates cell center location. The function (r) and Q are the Minmod flux limiter and the dissipation function, respectively, and they are defined as follows:
50 212112121123221,)],,1min(,0max[)( iiiiiiiiiiffffrffffrrr (3-57) bifbbifbbQQii,,21)(22121 . (3-58) The parameter is used to regulate the numerical viscosity; in this study it is assigned as zero, which results in the clipping of the fluxes. Local characteristic speed b is assigned according to the value of the cell face velocity. 3.4.1 Pressure-Based Algorithm for Steady-State Computations of Cavitating Flows The pressure-based algorithm, adopted for steady-state computations, follows the spirit of the well-established SIMPLE algorithm (Patankar and Spalding 1972, Patankar 1980), with substantial extension to treat issues associated with curvilinear coordinates and multiblock interface. Basically, the momentum equations are discretized as uPPdPnbunbPuPbPVuAuA )( , (3-59) where and are the coefficients of the cell center and neighboring nodes, respectively, due to contributions from convection and diffusion terms. V uPA unbA P and uPb represent the volume of the cell and the source term, respectively. Note that the d operator is the discrete form of the gradient operator. For the present study, there is no source term in the momentum equations; hence, it does not appear in the following formulation. The above form can also be written in a more convenient way (Moukalled and Darwish 2000) PdPPPPuu)(][ DH , (3-60)
51 where D P is written as wPPvPPuPPPAVAVAV/000/000/D . (3-61) For collocated grids, D P reduces to a scalar since. H is a linear operator resulting from the discretization of convection and diffusion terms. wPvPuPAAA The solution procedure is based on the predictor-corrector approach, where the discretized momentum equations is cast as PndPPPPDuu)(][1** H , (3-62) indicating that the velocity field at any given location is updated based on the existing values of the neighboring velocity and pressure. The superscript (n-1) stands for the values from the previous iteration. Substraction of Eq. (3-62) from Eq. (3-60) leads to the following correction term PdPPPPDuu)(][ H . (3-63) The steady-state continuity equation is written in discretized form as 0][ cfSnu . (3-64) In order to derive corrections for the pressure field, the continuity equation, Eq. (3-64), is converted to a pressure correction equation by substituting the corrected velocities: PPcfPdUSnPD])([* , (3-65) where cfSnuU . (3-66) Clearly, the above pressure correction equation has a diffusive nature. Note that the ][u H corrections are neglected to have a simpler form of the pressure-correction
52 equation. A variety of pressure-based algorithms have been proposed in years following the SIMPLE algorithm (Patankar and Spalding 1972). The differences among these algorithms are due to different approximations related to the ][u Hu term. Moukalled and Darwish (2000) have reviewed the various pressure-based algorithms and present them in a unified formulation. * PDPC In the pressure-based algorithm, the pressure-correction equation has been revised to achieve successful solutions for highly compressible flows (Shyy and Braaten 1988, Karki and Patankar 1989, Shyy 1994). This formulation will be described in the context of noncavitating flows with compressibility effects to motivate the present cavitating flow method. For highly compressible flows, density needs to be corrected to account for the strong pressure-density dependency. For such a formulation the flux terms in the continuity equation is uuuuuu *****))(( , (3-67) where u term is the mass flux entering the control volume. Starred variables represent the predicted value and primed variables represent the correction terms. The inclusion of the above equation, along with a relation that couples the density to pressure, leads to the following pressure correction equation: PPPC (3-68) PcfdPPPPcfdSnPUUPCSnPD])()([** . (3-69) By comparing Eqs. (3-65) and (3-69), one can see that the characteristic of the pressure-correction equation is altered from a pure diffusive nature to a mixed convective-diffusive nature in regions where density is a function of pressure. As discussed in Shyy (1994), the relative importance of the first and second terms in Eq. (3-69) depends on the
53 local Mach number; for low Mach number flows, only the first term prevails, while for high Mach number the second term becomes important. The Mach number dependency can be shown through the equation of state. The fourth term is a nonlinear second-order correction term and it can either be neglected or included in the source term to stabilize the computation in early iterations. In the cavitation models described earlier, a convection equation with pressure dependent source terms is solved to determine the density field. Because of this coupling between pressure and density, the pressure-correction equation needs to be reformulated even though the Mach number effect is not explicitly addressed in the model. Once the cavitation model is implemented into a pressure-based algorithm, the pressure-correction equation exhibits a convective-diffusive nature in cavitating regions and purely diffusive nature in the liquid phase. In the present algorithm, the following relation between density correction and pressure correction is introduced to establish the pressure-density coupling PPlPPC )1( , (3-70) where C is an arbitrary constant. It should be emphasized that the choice of this constant does not affect the final converged solution because of the nature of the pressure-correction equation. Different values of C simply lead to different paths for reaching the converged solution. However, the convective-diffusive nature of the pressure correction equation is directly affected by the choice. It can easily be shown that the ratio between the convective strength, the second term in Eq. (3-69), and diffusive strength, first term in Eq. (3-69), is directly related to C(1L ). In the cavity region, the liquid mass fraction decreases, and the pressure-correction equation is of a clear convective-diffusive nature.
54 On the other hand, in the liquid region, the pressure-correction equation returns to a purely diffusive type. Furthermore, it is found that a very large value for C can destabilize the computation in early stage of the iteration process. For this reason, it is suggested that C=O(1) be used. In the present computations, C=4 is adopted. It should be noted that the above discussion is based on normalized variables, with density of the liquid assigned to be 1. The present pressure-velocity-density coupling scheme is along a path responsive to the cavitation dynamics. It is found that the proposed scheme mimics the P variation adequately, especially near L =0, where the variation is very steep. From the convergence rate point of view, the scheme performs satisfactorily compared to the non-cavitating flow computations. Due to the convective-diffusive nature of the pressure-correction equation for cavitating flows, the coefficient matrix is noticeably asymmetric. Hence, iterative matrix solvers, designed specially for the fast efficient solution of the symmetric problems need new insight into the precondition treatment. In the present study, the conventional relaxation technique is employed. Next, the practical implementation of the pressure-velocity-density coupling for cavitating flows is presented. For simplicity, the following one-dimensional case is considered to illustrate the nature of the revised pressure correction equation. Note that a first order upwind scheme is used for the second term in Eq. (3-69): iiiilinciiiilinciiibPuCaPuCaPa1211112111]0,max[)1(]0,max[)1( (3-71)
55 ]0,max[)1(]0,max[)1(211211iilinciiilinciiuCauCaa . (3-72) In this equation, and a are the coefficients stemming from an incompressible formulation, b incia1 inci1 i is the source term, and L is the liquid volume fraction. Subscripts (i+1) and (i-1) stand for the neighboring grid nodes in the east and west direction, respectively. The above form is a combined incompressible-compressible formulation that preserves the incompressible nature in the liquid phase. In the cavitating region, it accounts for the pressure-density dependency in a nonlinear fashion, in accordance with the local value of L . This modification is key to a stable computation in which the uniform vapor pressure is recovered in the final converged solution. Another aspect is that, similar to compressible flow computations, the density at the cell face is upwinded both in the discretized momentum and pressure correction equations (Shyy 1994). The criterion for upwinding is based on the value of liquid volume fraction; that is, wherever L is less than 1.0, the cell-faced density value is estimated based on an upwinded formula. In regions of sharp density gradients, a single point upwinded extrapolation for density, instead of a two-point interpolation can significantly improve the convergence level. It should also be emphasized that Eq. (3-69) is not limited to the cavitation model employed in this study; it can easily be adopted for other cavitation models. Finally the major steps of the algorithm are summarized. 1. Compute the discretized momentum equation, Eq. (3-62). 2. Compute the convective-diffusive pressure-correction equation, Eq. (). 3-69 3. Correct the velocity and pressure field.
56 4. Compute the discretized scalar equations for L , k and , resulting from cavitation and turbulence models, using the corrected velocity and pressure field. 5. Compute the density field, Eq. (3-30), and turbulent viscosity field, Eq. (3-16). 6. Go to Step 1 and iterate until convergence. 3.4.2 Pressure-Based Algorithm for Time-Dependent Computations of Cavitating Flows In SIMPLE type of pressure-based methods, the equations are solved successively by employing iterations. The advantage gained by implicit differencing of the equations is coarsened during iterative process and time dependent computations can be highly expensive due to need for iterations. Issa (1985) have developed a pressure-based algorithm called the Pressure-implicit with Splitting of Operators (PISO) for the solution of unsteady flows. For steady flows this algorithm can be interpreted as an extension of the SIMPLE algorithm. PISO algorithm employs the splitting of operations in the solution of the implicitly discretized momentum and pressure or pressure-correction equations. This makes the solution procedure sequential in time domain and enables the exact solution of discretized equations at each time step with a formal order of accuracy. Operator splitting method eliminates the need for under relaxation adopted in SIMPLE type of algorithms. Issa (1985) has presented the compressible formulation of the algorithm. In this formulation the pressure-density coupling is introduced only through the time dependent term of the continuity equation. Bressloff (2001) has extended the PISO algorithm for high-speed flows by adopting the pressure-density coupling procedure in all-speed SIMPLE type of methods.
57 In the case of cavitating flow computations, the typical relaxation factors used in the iterative solution process are smaller compared to the ones used in computations of single-phase flows, and in addition smaller time steps are needed to study the cavitation dynamics. Hence, time-dependent computations with an iterative SIMPLE type pressure-based method become impractical. To eliminate this difficulty, the PISO algorithm, formulated in Thakur and Wright (2002) for curvilinear coordinates with all-speed treatment, is modified and extended for unsteady cavitating flow computations. The formulation of the present algorithm shares similar features of the formulation presented in Moukalled and Darwish (2000). In the predictor step the discretized momentum equations are solved implicitly using the old time pressure to obtain an intermediate velocity field. A backward Euler scheme is used for the discretization of the time derivative term. tuPDuunPndPPP11**)()(][ H (3-73) The intermediate velocity field does not satisfy continuity and needs to be corrected using the continuity equation as a constraint. In the first corrector step, a new velocity field,u, and a new pressure field, ** * P , are sought. The discretized momentum equation at this step is written as tuPDuunPPdPPP1****)()(][ H . (3-74) 3-74 Eq. (3-73) is substracted from Eq. (), leading to velocity correction terms PdPPPPDuu)(*** . (3-75)
58 If the pressure field depends on the density field, such as in high-speed flows or in cavitating flows, density field needs to be corrected and the following is written, as previously discussed in the previous section PnPP1* PPPC . (3-76) The discretized continuity equation written for the new velocity field and density field reads the following 0][***1*PcfPnPPSnuVt . (3-77) Inclusion of the corrected velocity field, Eq. () and density field, Eq. (3-76), transforms the above discretized continuity equation into a pressure-correction equation to be solved in the first corrector step 3-75 PnPPcfdnPPUUPCSnPDVtPC])([*1*1 . (3-78) When Eq. (3-73) is substracted from Eq. (3-74), the terms containing operator H is lost. Hence a second corrector step is needed to satisfy the mass conservation. In the second corrector step a new velocity field, ***Pu , and a new pressure field, ** P , are sought. The discretized momentum equation at this step is written as tuPDuunPPdPPP1*******)()(][ H . (3-79) Eq. (3-74) is substracted from Eq. (3-79), leading to the following correction terms PdPPPPPDuuuu)(][******** H . (3-80) The density field needs to be corrected at this step also as follows: PPnPPPP1*** . (3-81)
59 Inclusion of the corrected velocity field, Eq. () and density field, Eq. (3-81), transforms the discretized continuity equation into a pressure-correction equation to be solved in the second corrector step 3-80 PcfPPPcfdPPSnuuUUPCSnPDVtPC]][)(*[*********H . (3-82) Issa (1985) has shown that addition of each corrector step increases the accuracy by order of one in time. In this study, it is found that, instead of adding an additional corrector step, repeating the first corrector and second corrector step one more time gives the correct time-dependent behavior, such as the vortex shedding behind a cylinder. Additionally, since the splitting error depends on time step size, smaller time steps are needed for time accurate computations. Oliveria and Issa (2001) have discussed the coupling of temperature equation to the PISO algorithm for buoyancy-driven flows. Their study has shown that the PISO algorithm is amenable to different arrangements in the operator splitting procedure. An improved method is developed to handle the strong coupling of velocity and temperature in buoyancy-driven flows. For turbulent cavitating flow computations, special attention is also needed to couple the cavitation model and the turbulence model equations in the operator splitting procedure. The sequence of operations in PISO for turbulent cavitating flow computations is summarized below. 1. Predictor Step: Solve the dicretized momentum equation, Eq. (3-73), implicitly. 2. Compute needed for pressure-density coupling. C 3. First Corrector Step: Solve the convective-diffusive pressure-correction equation, Eq. (3-78), implicitly.
60 4. Correct the velocity and pressure field. 5. Solve implicitly the discretized scalar equations for L , k and , resulting from cavitation and turbulence models, using the corrected velocity and pressure field. 6. Compute the density field, Eq. (3-30), and turbulent viscosity field, Eq. (3-16). 7. Compute needed for pressure-density coupling. C 8. Second Corrector Step: Solve the convective-diffusive pressure correction equation, Eq. (3-82), implicitly. 9. Correct the velocity and pressure field. 10. Solve explicitly the discretized scalar equations for L , k and , resulting from cavitation and turbulence models, using the corrected velocity and pressure field (no need for a matrix solver). 11. Compute the density field, Eq. (3-30), and turbulent viscosity field, Eq. (3-16). 12. Go to Step 2 and repeat the steps one more time to enhance the coupling. 13. Proceed to the next time step. As explained in the previous section, the density at the cell face is upwinded based on a single point extrapolation, both in the discretized momentum and pressure-correction equations, to enhance mass and momentum conservation in regions of sharp density gradients. 3.4.3 The Issue of Speed of Sound Definition In Section 3.4.1, the importance of pressure-density coupling for cavitating flow computations have been discussed. The discussion still holds for the present algorithm. For steady-state computations, the pressure-density coupling scheme only affects the convergence path. The final solution is independent of the choice because of the nature of pressure-correction. However, for unsteady computations the choice of pressure-density coupling becomes critical, since it relates the propagation of information from the
61 cavitating region to the rest of the domain. The pressure-density coupling term C is a measure of the isentropic speed of sound. The relation is given as follows: PC 21cPCS , (3-83) where c is the speed of sound. In high-speed flows, the exact form of the speed of sound can be computed easily from the equation of state. However, in cavitating flows, computation of the speed of sound is not an easy task. Each transport equation-based cavitation model defines a different speed of sound as a result of a more complex functional relationship. In literature, there have been theoretical studies on defining the speed of sound in multiphase flows (Wallis 1969). One dimensional assumption and certain limitations are typical in these studies. Obviously, these definitions do not necessarily represent the actual speed of sound imposed by the cavitation models of interest in this study. They can only be an approximation. On the other hand, the fundamental definition of speed of sound as given in Eq. (3-83) could be useful, should the path to compute the partial derivative is known. From these arguments, it is clear that the computation of the speed of sound in cavitating flows is an open question. In this study, two different definitions for the speed of sound is investigated to be utilized in the pressure-density coupling scheme. The first one is the previous pressure-density coupling scheme used for steady-state computations as given below: )1(12LSCcPC . (3-84) The above form is referred to as SoS-1 in the rest of the study.
62 The second one is based on an approximation made to the fundamental definition of speed of sound. It is assumed that the path to compute the partial derivative P is the mean flow direction. This direction is preferred solely because the pressure variations are dominant in the mean flow direction (). This second definition is referred to as SoS-2 in the rest of the study and it is given as follows: 1111iiiiSPPPPC . (3-85) The partial derivative is computed based on central differencing of the neighboring nodes. The absolute value function is introduced to make sure a positive value is computed. 3.4.4 Boundary Conditions Velocity components, volume fractions and turbulence quantities are specified at the inlet boundary. At the outlet, a zero gradient condition is imposed for pressure correction, velocity and volume fraction; in addition, the velocity components are corrected to satisfy the global mass conservation condition with the pressure being fixed at a single upstream node. The pressure value at this node is used to compute the value of vapor pressure for the corresponding cavitation number. At walls, pressure, volume fractions and turbulence quantities are extrapolated along with no-slip conditions for velocity. For unsteady computations, a zero gradient condition for pressure correction causes difficulties because the global mass conservation condition involves the time derivative terms of density; hence, the pressure is specified on the outlet boundary for unsteady computations (Patankar 1980). The outlet boundary is placed far away from the cavitating region of interest. This type of outlet boundary condition does not require a
63 global mass conservation condition since the pressure correction is exactly zero because of the specified pressure. Boundary conditions for other variables are the same as in steady-state computations.
CHAPTER 4 RESULTS OF STEADY-STATE COMPUTATIONS 4.1 Validation of the Pressure-Based Algorithm for Cavitating Flows The pressure-based algorithm for steady-state computations of cavitating flows, presented in the previous chapter, is applied to two different axisymmetric geometries to investigate numerical and physical issues involved in cavitation modeling. Both geometries considered have a cylindrical after-body and they will be referred to according to the shape of their headforms, specifically, hemispherical or blunt. Model-2 (Kunz et al. 2000) is adopted for cavitation modeling. The results include steady-state computations of noncavitating and cavitating flows at a Reynolds number of 1.36x10 5 . Results are compared with experimental results of Rouse and McNown (1948), in which pressure coefficients are reported. Since the experimental information does not report time dependency, the steady-state model is adopted in the present computations. From the physical point of view, the steady-state assumption is sensible for sheet cavitation, which has a quasi-steady behavior, with most of the unsteadiness localized in the rear closure region (Knapp 1970, Brennen 1995). However, for the blunt object, cavitation may not develop in the form of sheet cavitation. Nevertheless, the steady-state formulation is employed for both geometries. The results are organized in two parts based on the geometry considered. First, the simulations of hemispherical object are presented and issues of grid refinement, convection schemes, and sensitivity to cavitation model parameters are discussed. Next, 64
65 computations of the blunt object are presented in which the turbulence modeling issue in the context of the kmodel is investigated. Information related to velocity, pressure and cavity characteristics are discussed. 4.1.1 Flow over a Hemispherical Object A grid refinement study has been performed to assess the accuracy and sensitivity of the predictions for flows around a hemispherical object. Two grid systems (Grid-A with 119x65 nodes; Grid-B with 231x161 nodes) are generated. Compared to Grid-A, Grid-B has essentially twice the spatial resolution in each direction surrounding the cavity. The computational domain, corresponding to Grid-A, and the imposed boundary conditions are shown in Figure 4-1. In both computations, the same modeling parameters, given in are used. shows the pressure coefficient, flow pattern and density distribution along the body obtained at a cavitation number of 0.40. As can be seen from the plot, the impact of grid resolution on the pressure distribution is not significant. Also no sharpening of the interface is observed due to grid refinement. However, the density profiles indicate a shorter cavity length for Grid-B. The reason for this shortening is primarily due to the reentrant jet in the rear closure region. By comparing the streamlines and cavity profiles, given also in , one can see that the reentrant jet in the closure region is clearer with Grid-B, and has a higher strength than with Grid-A. Furthermore, a lower density is observed inside the cavity with Grid-B. It should be emphasized that the pressure profile is not sensitive to the range of grid resolution investigated. Although the flow structure is better resolved with Grid-B, Grid-A is utilized for the rest of the study because of the negligible impact on pressure distribution. Figure 4-2Figure 4-2Figure 4-2
66 This decision is made based on the fact that only pressure distributions are available in the experimental study by Rouse and McNown (1948). Figure 4-1. The computational domain and imposed boundary conditions. 0 0.5 1 1.5 2 2.5 3 3.5 4 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/dCp, densityHemispherical, =0.40Re=1.36x105 Grid-A (119x65) Grid-B (231x161)Exp. data Grid-A Grid-B density 0 0.5 1 1.5 2 2.5 3 3.5 4 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/dCp, densityHemispherical, =0.40Re=1.36x105 Grid-A (119x65) Grid-B (231x161)Exp. data Grid-A Grid-B density Grid-A (119x65) minl=47 Grid-B (231x161) minl=114 Figure 4-2. Grid refinement study for hemispherical object with =0.40 (C dest =9x10 5 , C prod =3x10 4 , L / V =1000).
67 0 500 1000 1500 2000 2500 3000 -6 -5 -4 -3 -2 -1 0 1 iteration numberlog(residual) P'lno upwindingupwinding Hemispherical, =0.40,Re=1.36x105lP' 0 500 1000 1500 2000 2500 3000 -6 -5 -4 -3 -2 -1 0 1 iteration numberlog(residual) P'lno upwindingupwinding Hemispherical, =0.40,Re=1.36x105lP' Figure 4-3. Improvement in convergence level due to density upwinding. Hemispherical object at =0.40. Figure 4-3 As explained previously in the numerical methodology section, the density at the cell face is upwinded both in discretized momentum and pressure-correction equations. In regions of sharp density gradients, a single point upwinded extrapolation for density can significantly improve the convergence level as demonstrated in . The residuals resulting from momentum, mass continuity, and volume fraction transport equations can approach much lower levels than the pressure-based algorithm without an upwinded interpolation for density. It is noted that the residuals are defined as the absolute values of the imbalance of the individual equations summed over the entire number of computational cells, normalized by the total flux of the given variable, at the inlet of the computational domain. Numerical diffusion is a major concern in complex flow computations. In this study, the effect of first and second order accurate convection schemes on solution accuracy is
68 evaluated. In Figure 4-4, the pressure coefficient and density profiles, corresponding to a cavitation number of 0.40, are shown. Interestingly, both convection schemes produce almost identical results for density and pressure. The density profile is very sharp at the closure region even in the case of first order upwind scheme (FOU). This phenomenon is mainly an outcome of the cavitation model. Basically the variable density field is generated through source terms for destruction and production of liquid phase. Any smearing is automatically eliminated through activation of the production source term and as the iteration proceeds a balance between production and destruction terms is achieved resulting in a cavity profile with uniform pressure and density field inside. No significant change is observed in flow structure and cavity profiles. Also the second order accurate CVS scheme does not produce any oscillations around the sharp discontinuity at the closure region. Even though the first order solution appears satisfactory in this case, to ensure that satisfactory results are obtained in other cases, the second order accurate CVS scheme is adopted for convection terms throughout the study. The sensitivity of the solutions to model parameters is studied in . It can be seen that even increasing these parameter by an order of magnitude has little effect on the pressure coefficient predictions. However, the computed density ratio is noticeably different between these model parameters. Clearly, the computed density ratios can be controlled through adjustment of the model parameters to yield very different solutions while pressure predictions remain little unaffected. Figure 4-5
69 0 0.5 1 1.5 2 2.5 3 3.5 4 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/dCp, densityHemispherical, =0.40, Re=1.36x105 Cp, FOU Cp, CVS , FOU , CVS Cp, Exp. data 0 0.5 1 1.5 2 2.5 3 3.5 4 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/dCp, densityHemispherical, =0.40, Re=1.36x105 Cp, FOU Cp, CVS , FOU , CVS Cp, Exp. data Figure 4-4. Effect of convection schemes on pressure coefficient and density distribution for the hemispherical object at =0.40. (C dest =9x10 5 , C prod =3x10 4 , L / V =1000). Grid-A (119x65) is used. 0 0.5 1 1.5 2 2.5 3 3.5 4 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/dCp, densityHemispherical, =0.40,Re=1.36x105density Cdest=2x106, Cprod=1.3x105Exp. data Cdest=9x105, Cprod=3x104 0 0.5 1 1.5 2 2.5 3 3.5 4 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/dCp, densityHemispherical, =0.40,Re=1.36x105density Cdest=2x106, Cprod=1.3x105Exp. data Cdest=9x105, Cprod=3x104 Cdest=2x106, Cprod=1.3x105 minl=739 Cdest=9x105, Cprod=3x104 minl=47 Figure 4-5. Sensitivity of modeling parameters for the hemispherical object at =0.40 ( L / V =1000). Grid-A (119x65) is used.
70 0 1 2 3 4 5 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/dCpHemispherical, Re=1.36x105 noncavitating =0.40 =0.30 Exp.noncavitatingExp. =0.40 Exp. =0.30 0 1 2 3 4 5 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/dCpHemispherical, Re=1.36x105 noncavitating =0.40 =0.30 Exp.noncavitatingExp. =0.40 Exp. =0.30 =0.30 minl=653 =0.40 minl=47 Figure 4-6. Comparison of pressure coefficient distributions for hemispherical object under noncavitating and cavitating conditions (C dest =9x10 5 , C prod =3x10 4 , L / V =1000). Figure 4-6Figure 4-6 demonstrates the predictive capability of the model at cavitation numbers of 0.40 and 0.30 through comparison with experimental data (Rouse and McNown 1948). Identical model parameters are adopted for both cavitation numbers. The pressure distribution corresponding to the noncavitating condition is also plotted for comparison. The present numerical algorithm performs well for both cavitating and noncavitating conditions. The corresponding cavity profiles, streamlines and computed density ratios are also presented in . The computed cavity profiles are in the form of pinched pockets with reentrant jets in the closure region. With a lower cavitation number (=0.30), the cavity, as expected, becomes larger than that with a cavitation number (=0.40). The reentrant jet is also stronger suggesting that at lower cavitation numbers the reentrant jet can easily perturb the cavity, possibly leading to shedding of bubbles. The computed density ratio is higher for =0.30 because, similar to what is observed in grid refinement, the source terms are effective on more grid points. The cavity
71 detachment point remains fixed in both of the simulations, which is also in agreement with experimental data. Figure 4-7. Vorticity generation at the closure region due to baroclinic torque. Gopalan and Katz (2000) have recently shown that vorticity production occurs at the closure region of sheet cavities due to baroclinic torque. As shown in Figure 4-7, the present computations do indicate production of vorticity at the closure region for both cavitation numbers considered, which is consistent with the findings of Gopalan and Katz (2000) and, hence, verifies the success of the present numerical methodology. Furthermore, there is no additional production of vorticity at the front part of the cavity indicating that the density and the pressure fields are properly computed so that their gradients align parallel to each other without causing any baroclinic vorticity generation. The dynamics of the phase change process is shown in Figure 4-8. The evaporation process, driven by the pressure difference, is localized at the leading part of the cavity generating the vapor phase. On the other hand, condensation is concentrated along the interface. As a result, these two counteracting processes generate a vapor pocket, which has an almost uniform density (as shown in Figure 4-7) and pressure field inside.
72 =0.40=0.30Evaporation contoursEvaporation contoursCondensation contoursCondensation contours =0.40=0.30Evaporation contoursEvaporation contoursCondensation contoursCondensation contours Figure 4-8. Dynamics of the phase change. 4.1.2 Flow over a Blunt Object Computations have also been conducted to predict cavitation over blunt cylindrical objects. The computational domain and the imposed boundary conditions are shown in . For this particular object, Rouse and McNown (1948) have provided the pressure coefficient distribution along the body, and Katz (1984) has reported flow reattachment locations. In view of the potential deficiencies associated with the original kturbulence model (Jones and Launder 1972) for complex flows, a nonequilibrium version of the kmodel (Shyy et al. 1997) is studied in addition to the original kmodel. To investigate the turbulence modeling issue, the noncavitating condition is considered first and then the study is extended for cavitating conditions. Figure 4-10 shows the pressure coefficient profile for the noncavitating condition. The original kmodel fails to match the experimental data in the vicinity of the sharp corner, where large strains and streamline curvatures are expected to occur. On the other hand, the nonequilibrium kmodel performs better to capture the pressure coefficient distribution. Both models produce comparable solutions far downstream. It should be pointed out that Figure 4-9
73 at far downstream, pressure distributions are slightly different compared to the experimental data. This is apparently because of the difference between the computational domain and the experimental test section. Figure 4-9. The computational domain and imposed boundary conditions. Figure 4-11 shows the normalized u-velocity distribution along the body that is extracted from the first computational nodes away from the solid boundary. The point where velocity changes sign is the reattachment point. One can see that there are considerable differences in predictions of reattachment locations and associated velocity profiles between two turbulence models. It is observed that the reattachment point location, predicted with the nonequilibrium kmodel, is consistent with the experimental data (Katz 1984) even though the Reynolds number range is not the same. It should be noted that Katz (1984) has not observed a strong correlation between
74 reattachment location and Reynolds number for this range. Overall, the nonequilibrium kmodel performs better, compared to the original kmodel, in capturing the pressure coefficient distribution and the extent of separation zones. Based on this observation, the nonequilibrium kmodel is included in the study of cavitating flows over blunt objects. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -1.5 -1 -0.5 0 0.5 1 1.5 s/dCpBlunt, Single-phase, Re=1.36x105 Standard k-Noneq. k-Exp. data 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -1.5 -1 -0.5 0 0.5 1 1.5 s/dCpBlunt, Single-phase, Re=1.36x105 Standard k-Noneq. k-Exp. data Ori g inal Figure 4-10. Impact of turbulence modeling on predictions of pressure coefficient distribution for noncavitating conditions. Figure 4-12Figure 4-12 shows the pressure coefficients and density profiles along the blunt object at a cavitation number of 0.50. Cavity profiles and streamlines are also included in . Computed pressure profiles qualitatively follow the experimentally observed trend. Unlike noncavitating flow computations, both turbulence models considered produce similar results. illustrates the normalized u-velocity profile along the surface. The difference in predicted reattachment location is also less significant between the two turbulence models. Furthermore, both solutions exhibit rather modest density ratios across the cavity interface, suggesting that inside the cavity the flow exhibits liquid-vapor bubbly structures. In our view, the discrepancy between computational and experimental Figure 4-13
75 results may be caused by our steady-state computations. Stinebring et al. (2001) employ high speed visualization and report that cavities formed around blunt objects are highly unsteady and bubbling phenomenon is typically observed. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 s/dU/UoBlunt, Single-phase, Re=1.36x105 Standard k-Noneq. kRange of reattachment points for Re~2x105-1x106Katz (1984) Present computationrecirculationzonestandard kPresent computationrecirculationzonenonequilibriumk0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 s/dU/UoBlunt, Single-phase, Re=1.36x105 Standard k-Noneq. kRange of reattachment points for Re~2x105-1x106Katz (1984) Present computationrecirculationzonestandard kPresent computationrecirculationzonenonequilibriumkOri g inal Figure 4-11. Effect of turbulence model on separation zones for noncavitating conditions. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -1.5 -1 -0.5 0 0.5 1 1.5 s/dCp, densityBlunt, =0.50, Re=1.36x105 Cp, Standard k-Cp,Noneq. k-, Standard k-, Noneq. k-Exp. data 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -1.5 -1 -0.5 0 0.5 1 1.5 s/dCp, densityBlunt, =0.50, Re=1.36x105 Cp, Standard k-Cp,Noneq. k-, Standard k-, Noneq. k-Exp. data Nonequilibrium kminl=2 Standard kminl=6 Ori g inal Ori g inal Ori g inal Figure 4-12. Impact of turbulence modeling on predictions of pressure coefficient and density distribution (=0.50, C dest =9x10 3 , C prod =6x10 3 , L / V =1000).
76 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 s/dU/UoBlunt, =0.50, Re=1.36x105 Standard k-Noneq. kReattachmentpoints 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 s/dU/UoBlunt, =0.50, Re=1.36x105 Standard k-Noneq. kReattachmentpoints Ori g inal Figure 4-13. Comparison of u-velocity along the surface obtained from different turbulence models. (=0.50, C dest =9x10 3 , C prod =6x10 3 , L / V =1000). In Figure 4-14, the turbulent kinetic energy k and turbulent dissipation rate , resulting from the original and nonequilibrium kturbulence models, are compared. In both cases, the turbulent variables are of higher levels following the shear layer surrounding the cavity; however, the detailed distributions are different. The nonequilibrium model yields more substantial presence in both k and . Together, they are responsible for a somewhat longer recirculating zone, as shown in Figure 4-13. Although the predictions of pressure are less sensitive to turbulence modeling, other quantities such as wall shear stress and velocity profiles can be more affected by it. Due to lack of experimental guidance, detailed results from the computational study are not presented.
77 Original k k turbulent kinetic energy N onequilibriu m k-, k turbulent kinetic energy Original k-turbulent dissipation rate N onequilibriu m k-, -turbulent dissipation rate k k turbulent kinetic energy k k turbulent kinetic energy Figure 4-14. Comparison of turbulence quantities obtained from different turbulence models. (=0.50, C dest =9x10 3 , C prod =6x10 3 , L / V =1000). 4.2 Evaluations of Transport Equation-Based Cavitation Models In this section, empirical transport equation-based cavitation models, summarized in Section 3.3.1 as Model-1 (Singhal et al. 1997, Merkle et al. 1998), Model-2 (Kunz et al. 2000), Model-3 (Singhal et al. 2001), and the new interfacial dynamics model, developed in Section 3.3.2, are assessed using the unified incompressible-compressible pressure-based algorithm. The new interfacial dynamics cavitation model is further studied to investigate the flow structure of attached stable cavities. The original kmodel is adopted for turbulence closure (Jones and Launder 1972). Three flow configurations have been considered, namely, (i) an axisymmetric cylindrical object with a hemispherical forehead (referred to as hemispherical object) at a Reynolds number of 1.36x10 5 and a cavitation number of 0.30, (ii) the NACA66MOD hydrofoil at an angle of attack of 4 with Reynolds number of 2.0x10 6 , and (iii) the
78 convergent-divergent nozzle at a Reynolds number of 2.0x10 6. It is experimentally observed that sheet (attached) cavitation occurs for all three of the geometries under given conditions and time averaged experimental data of pressure distribution along the surface is available for hemispherical object and the NACA66MOD hydrofoil (Rouse and McNown 1948, Shen and Dimotakis 1989). Stutz and Reboud (1997a) have detailed the two-phase flow structure of stable sheet cavitation in the convergent-divergent nozzle and have provided time averaged velocity and vapor volume fraction profiles within the attached cavity. 4.2.1 Flow over a Hemispherical Object Figure 4-15 compares the performance of the cavitation models for the hemispherical object at a cavitation number of 0.30. The computational domain and the imposed boundary conditions have been shown in Figure 4-1 previously. All three models match the experimental data of pressure coefficient distribution satisfactorily. Differences in the performance are more noticeable in the closure region, where the vapor phase condenses. shows the corresponding density distribution along the surface. As seen from density plots, the liquid phase first expands and vapor phase appears uniformly inside the cavity, then the vapor phase compresses, in a shock like fashion, back to the liquid phase. The differences among the three models in density profiles are significant. This implies that each model generates different compressibility characteristics, although they produce very similar steady-state pressure distributions. This issue has important implications in time dependent problems; which model produces the correct compressibility is an open question and needs further investigation. Figure 4-16
79 0 0.5 1 1.5 2 2.5 3 3.5 4 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/DCpHemispherical, Re=1.36x105, =0.30 Model-1 Model-2 Model-3 Exp. data 0 0.5 1 1.5 2 2.5 3 3.5 4 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/DCpHemispherical, Re=1.36x105, =0.30 Model-1 Model-2 Model-3 Exp. data Figure 4-15. Comparison of surface pressure distribution obtained from empirical cavitation models. Experimental data is from Rouse and McNown (1948). The performance of the new interfacial dynamics cavitation model is assessed through comparison with its empirical counter part, Model-1. Similar to what is observed in Figure 4-15, the prediction of pressure remains insensitive to the cavitation model choice. A noticeable difference exists in density distribution within the cavity. The density ratio is higher, close to 1000, for the case of the new interfacial dynamics model. Figure 4-18Figure 4-18 shows the distribution of density throughout the cavity and the spanwise vorticity distribution obtained from all four cavitation models. The interface is captured sharply with the new interfacial dynamics model, Model-2 and Model-3 compared to Model-1, especially in the downstream region of the cavity. As seen from the spanwise vorticity distribution, also given in , all models capture the additional generation of vorticity at the closure region. It is clear that vorticity generation in the
80 closure is not due to model choice but rather because of the fundamental modeling approach that is the transport equation modeling coupled to the unified pressure-based algorithm developed in this study. The model choice, among the transport equation based cavitation models, mainly affects the magnitude of the vorticity. 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.2 0.4 0.6 0.8 1 s/DdensityHemispherical, Re=1.36x105, =0.30 Model-1Model-2Model-3 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.2 0.4 0.6 0.8 1 s/DdensityHemispherical, Re=1.36x105, =0.30 Model-1Model-2Model-3 Figure 4-16. Comparison of surface density distribution obtained from empirical cavitation models.
81 0 0.5 1 1.5 2 2.5 3 3.5 4 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/DCp,densityHemispherical, Re=1.36x105, =0.30 Model-1, Cp Model-1, Exp. data 0 0.5 1 1.5 2 2.5 3 3.5 4 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/DCp,densityHemispherical, Re=1.36x105, =0.30 New Model, Model-1, Cp Model-1, Exp. data New Model, Cp 0 0.5 1 1.5 2 2.5 3 3.5 4 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/DCp,densityHemispherical, Re=1.36x105, =0.30 Model-1, Cp Model-1, Exp. data 0 0.5 1 1.5 2 2.5 3 3.5 4 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 s/DCp,densityHemispherical, Re=1.36x105, =0.30 New Model, Model-1, Cp Model-1, Exp. data New Model, Cp Figure 4-17. Comparison of surface pressure and density distribution obtained from empirical Model-1 and the new interfacial dynamics model. Experimental data is from Rouse and McNown (1948). Figure 4-18. Detailed view of the cavitating region. On the left is the density distribution; on the right is the spanwise vorticity distribution
82 4.2.2 Flow over a NACA66MOD Hydrofoil The NACA66MOD hydrofoil case is computed at a Reynolds number of 2x10 6 . The turbulent boundary layer is extremely thin at such high Reynolds numbers. Since the original kturbulence model along with the wall function is adopted, it is important to offer spatial resolutions consistent with the modeling requirement as previously studied for a single-phase flow over a flat plate in Section 3.2.3. This requires that the non-dimensional normal distance from the wall (y + ), a representation of the local Reynolds number, should be in the log-law region. Once a cavity occurs on the surface, the local Reynolds number decreases due to reduction in density and the first grid point away from the wall may not be positioned in the log layer. This issue has important implications in both accuracy and convergence. It is found that proper grid distribution is very important for satisfactory results and convergence of cavitating flow computations. However, there is a trade-off between positioning the first grid point away from the wall in the log layer verses having enough points to discretize the cavity, especially for high Reynolds number cases. Two grids have been tested. Model-2 is used in the computations. Grid-A is a three block domain with approximately 6.6x10 4 nodes, whereas Grid-B is a six block domain with approximately 9.8x10 4 nodes. Close up views of the same region in both grids are shown in Figure 4-19 to demonstrate the spatial resolution. The boundary of the cavitating region is also highlighted on the grids to give an idea of how many grid points have been used to discretize the vapor cavity. Note that the first grid point in Grid-B is positioned to satisfy the (y + ) requirement and the rest of the points are clustered close to it. However the grid generation software has not enabled further clustering of the grid points above the first point due to its own restrictions.
83 X Y -0.23 -0.225 -0.22 -0.215 -0.21 0.02 0.025 0.03 0.035 0.04 X Y -0.23 -0.225 -0.22 -0.215 -0.21 0.02 0.025 0.03 0.035 0.04 X Y -0.23 -0.225 -0.22 -0.215 -0.21 0.02 0.025 0.03 0.035 0.04 X Y -0.23 -0.225 -0.22 -0.215 -0.21 0.02 0.025 0.03 0.035 0.04 Figure 4-19. Visual compassion of the grid distributions in the cavitating region 0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 1.5 X/C-CpNACA66MOD, Re=2.0x106, =0.91 Grid-A Grid-B Exp. data 0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 1.5 X/C-CpNACA66MOD, Re=2.0x106, =0.91 Grid-A Grid-B Exp. data Figure 4-20. Effect of grid distribution on predictions. Model-2 is used. Experimental data is from Shen and Dimotakis (1989).
84 0 0.2 0.4 0.6 0.8 1 100 101 102 103 104 X/CY+NACA66MOD, Re=2.0x106, =0.91 Grid-AGrid-B y+=11.63 y+=400.0 0 0.2 0.4 0.6 0.8 1 100 101 102 103 104 X/CY+NACA66MOD, Re=2.0x106, =0.91 Grid-AGrid-B Grid-AGrid-B y+=11.63 y+=400.0 Figure 4-21. Distribution of the y + values of the first grid points away from the wall along the hydrofoil surface As shown in Figure 4-20 the predictions with Grid-A are very poor compared to the case with Grid-B. A very short cavity has been captured with Grid-A as a result of poor resolution. The vapor cavity is so thin that Grid-A can only accommodate two grid points in this region. In Figure 4-21 the y + distribution of the first grid points away from the wall is plotted along the surface. The law of the wall is also highlighted on this plot. As can be observed in Figure 4-19, special attention has been given during grid generation to place the points of the Grid-B in the log layer. A higher resolution grid is also utilized and it is found that if the y + values of the cavitating region are in the viscous sublayer, the computation is not stable and convergent. This may be because in the wall function formulation a linear velocity profile is imposed in the viscous sublayer and such a profile may not be suitable to represent the phase change dynamics. It should be emphasized that
85 the original kturbulence model is based on the equilibrium assumption, and the rate of turbulence production and dissipation are balanced in the log-layer; hence, it is advisable to place the near wall region in the log layer. If higher resolution is required a low Reynolds number two-equation turbulence model should be the choice. Nevertheless Grid-B is used for all other computations of NACA66MOD hydrofoil case. In Figure 4-22 and Figure 4-23, the performances of the cavitation models have been assessed for cavitation numbers of 0.91 and 0.84, respectively. For cavitation number of 0.84, the vaporous cavity covers about 60% of the hydrofoil surface. An overall flow visualization is presented in Figure 10. All three models produce satisfactory results at both cavitation numbers. Differences in predictions are more pronounced, similar to the hemispherical object case, at the closure region, which is due to the different compressibility characteristics imposed by the cavitation models. Model-2 and Model-3 do relatively a better job to capture the pressure distribution at the closure region as compared to the predictions of Model-1 In Figure 4-25, the result obtained from Model-1 is compared with results of the published studies in literature that utilize the same model also (Singhal et a 1997, Ahuja et al. 2001). Differences between each study are the numerical methods and the computational grids that are adopted. The present results show a better agreement with the experimental data as compared to other computational studies.
86 Figure 4-22. Comparison of surface pressure distribution. The cavitation number is 0.91. Experimental data is from Shen and Dimotakis (1989).
87 Figure 4-23. Comparison of surface pressure distribution. The cavitation number is 0.84. Experimental data is from Shen and Dimotakis (1989). Figure 4-24. Detailed flow visualization. Cross sectional view shows pressure distribution.
88 0 0.2 0.4 0.6 0.8 1-0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 X/C-CpNACA66MOD, Re=2.0x106, =0.91 Present study Ahujaet al. (2001) Singhalet al. (1997)Exp. data 0 0.2 0.4 0.6 0.8 1-0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 X/C-CpNACA66MOD, Re=2.0x106, =0.91 Present study Ahujaet al. (2001) Singhalet al. (1997)Exp. data Figure 4-25. Comparison of the present results obtained from Model-1 with the results of published studies in literature that also utilize Model-1. 4.2.3 Convergent-Divergent Nozzle (Stable Sheet Cavitation) There exists many experimental studies on cavitation but few of such studies have been useful to validate numerical studies. Most studies are generally limited to external liquid flows and mainly provide the pressure measurements. As shown previously, the transport equation based cavitation models, coupled to the unified pressure-based algorithm, have been successful to capture the pressure distributions on test geometries. To further investigate the predictive capability of the computational methodology, the two-phase flow structure within the cavity, the impact of cavitation on turbulence and vorticity dynamics are studied. Stutz and Reboud (1997a) have conducted experiments to detail the flow structure within sheet cavitation. A convergent-divergent nozzle is built in
89 a venturi type test section. The test section is 520 mm long 44 mm wide and 50 mm high at the inlet. The throat height is 43.7 mm. The angle of the convergent part is 4.3, whereas the divergent section angle is 4.0. A sharp ridge is formed at the throat section to induce further pressure drop. The experiments have been carried for a cavitation number ranging between 0.60-0.75. The Reynolds number is ranging from 4.3x10 5 to 2.1x10 6 . The accuracy of the velocity measurements is 10%. The results have been presented for a mean cavity length of 80 mm. The LDA measurements have indicated a nearly 2-D flow. Based on this information, the cavitation number and the Reynolds number are chosen as 0.62 and 2.0x10 6 , respectively. The computational domain is 2-D with minimum resolution in the third direction (infinitely wide nozzle assumption). The computational domain and the imposed boundary conditions are shown in Figure 4-26. Figure 4-26. The computational domain and the imposed boundary conditions for the convergent-divergent nozzle. Figure 4-27 shows a comparison of the noncavitating and cavitating conditions in the nozzle. The turbulent kinetic energy and the pressure fields for both cavitating and noncavitating conditions are shown. It is clear that cavitation not only alters the local flow structure but also alters the overall flow structure in the nozzle. The highest level of turbulence production is just downstream of the cavity. This finding is in good agreement with the experimental results of Gopalan and Katz (2000) who have also studied sheet
90 cavitation in a nozzle. It should be noted that the geometry and flow conditions in their study are different than the ones in this study. The pressure field in the throat section is noticeably affected due to the presence of sheet cavitation, also shown in Figure 4-27. As can be deduced from the contour plots, the pressure gradient field under cavitating conditions is more complex compared to the noncavitating conditions. Strong adverse pressure gradients are responsible for the reverse flow in the closure region, which does not exist in noncavitating conditions. Figure 4-27. The impact of cavitation on the structure of internal flow through the nozzle. It is important to know the source of turbulence production in the closure region. In particular, this information is helpful to identify possible areas of improvement in computational modeling of cavitation. It is worthwhile to note that the present study utilizes the original kmodel for turbulence closure. This model has well-known deficiencies for complex flow structure (Wilcox 1993). Gopalan and Katz (2000) concluded that the origin of the turbulence production is the collapse of the cavities in the closure region, which results in vorticity generation. Obviously, the collapse of the cavities can not be well simulated within the present steady-state assumptions but the
91 spanwise vorticity field has indicated additional vorticity in the closure region, as shown previously in the hemispherical object case. Considering the fact that an adverse pressure gradient generates a recirculation zone behind sheet cavitation, the main question to be asked is; what are the effective terms in the vorticity transport equation responsible for the generation of spanwise vorticity in the closure region. The vorticity transport equation is written as (Lesieur 1993): 221)(Put (4-1) Figure 4-28 shows the contour plots of the effective terms of the vorticity transport equation along with the spanwise vorticity distribution. Clearly, the baroclinic term of the vorticity transport equation is important in the closure region of the cavity. Hence, it can be concluded that the additional turbulence production is driven by the baroclinic vorticity generation. Figure 4-29 shows a direct comparison of the computed cavity and velocity vectors with the experiment. The results are in good agreement with the experiment and indicate the thickening of the boundary layer in cavitating regions. The boundary layer follows the cavity boundary but slightly thicker than the cavity. The experiments show the reverse flow behind the cavity but it is not presented in this particular plot. A more direct comparison of the cavity boundary and velocity profiles within cavity is also performed. As seen in Figure 4-30, the cavity boundary and it is overall shape is well captured with the computations. Discrepancy exists in terms of capturing the extent of the reverse flow. The present results indicate a shorter one, which results in different velocity profiles. Later in Section 5.2.3 it will be shown that time-averaged simulations produce much better agreement with experimental data.
92 Figure 4-28. The effective terms of the vorticity transport equation for cavitating conditions.
93 Figure 4-29. The velocity field around sheet cavitation. The bottom plot is reproduced from the experimental study of Stutz and Reboud (1997a) with permission from the authors.
94 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 0.12 X/Lcav(%) Y/Lcav(%) 010 0 0.2 0.4 0.6 0.8 1 1.2 0 0.1 0.2 X/LcavY/Lcav Exp.data Present computation 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 0.12 X/Lcav(%) Y/Lcav(%) 010 0 0.2 0.4 0.6 0.8 1 1.2 0 0.1 0.2 X/LcavY/Lcav 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 0.12 X/Lcav(%) Y/Lcav(%) 010 0 0.2 0.4 0.6 0.8 1 1.2 0 0.1 0.2 X/LcavY/Lcav Exp.data Present computation 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 0.12 X/Lcav(%) Y/Lcav(%) 010 0 0.2 0.4 0.6 0.8 1 1.2 0 0.1 0.2 X/LcavY/Lcav Exp.data Present computation 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 0.12 X/Lcav(%) Y/Lcav(%) 010 0 0.2 0.4 0.6 0.8 1 1.2 0 0.1 0.2 X/LcavY/Lcav 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 0.12 X/Lcav(%) Y/Lcav(%) 010 0 0.2 0.4 0.6 0.8 1 1.2 0 0.1 0.2 X/LcavY/Lcav Exp.data Present computation Figure 4-30. Comparison of the cavity boundary and the velocity profiles within the cavity. Experimental data is from Stutz and Reboud (1997a).
95 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 X/Lcav(%)Y/LcavVapor-volume fraction profiles 0100V(%) Exp. data Present computation 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 X/Lcav(%)Y/LcavVapor-volume fraction profiles 0100V(%) Exp. data Present computation Figure 4-31. Vapor volume fraction profiles within the cavity. Experimental data is from Stutz and Reboud (1997a). Figure 4-31 presents the comparison of vapor volume fraction distribution. The experimental data indicate a mixture in the reverse flow whereas the computation captures a pure liquid reverse flow. The vapor content is much higher in the computational results. As will be shown in Section 5.2.3 time-averaged simulations produce much better agreement with experimental data.
CHAPTER 5 RESULTS OF TIME-DEPENDENT COMPUTATIONS 5.1 Single-Phase Flow over a Cylinder To validate the pressure-based algorithm for time-dependent computations, the unsteady problem of flow over a cylinder is studied. The Reynolds number investigated is 100, at which the flow is laminar and exhibits the well-known vortex shedding phenomenon. A perturbation is given to the initial flow field to quickly reach the periodic oscillatory flow pattern. The influence of the imposed boundary conditions on the time-dependent behavior is studied. The computational domain and the imposed boundary conditions are shown in Figure 5-1. Williamson (1989) has performed experiments with an open domain. On the left boundary, the steady-state inlet condition with constant horizontal velocity is adopted. On the right, the outlet condition employing zero gradient extrapolation of the velocity and specified pressure is invoked. Several conditions on the top and bottom boundaries are tested. The inlet boundary condition imposes the velocity, no-slip boundary condition sets all the velocity components to zero, the slip boundary condition imposes a zero gradient for the velocity, and the outlet boundary condition specifies the pressure. presents the Strouhal number corresponding to each case. With the same boundary conditions, the computed and measured values of the Strouhal number are in excellent agreement. Table 5-1 96
97 Figure 5-1. Computational domain and imposed boundary conditions Table 5-1. Influence of top and bottom boundary conditions on the Strouhal number. Strouhal number (fD/U ) No-slip 0.1267 Slip 0.1658 Inlet 0.1750 Outlet 0.1640 Exp. (Williamson 1989) 0.1640 Figure 5-2 shows the periodic behavior of the cross-stream (v) component of the velocity obtained by probing at a single point downstream of the cylinder. The v component of the velocity alternates due to the vortex shedding behind the cylinder. Vortices form on the surface of the cylinder and shed downstream forming the well known von Karman vortex street as shown in Figure 5-3.
98 0.6 0.8Flow over a cylinder, Re=100 v-velocity 0.6 0.8Flow over a cylinder, Re=100 v-velocity 0 10 20 30 40 50 60 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 t/T v/U 0 10 20 30 40 50 60 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 t/T v/U Figure 5-2. Periodic oscillation of the cross-stream (v) component of velocity. Figure 5-3. Spanwise vorticity distribution. Snapshot of the vortex shedding and the von Karman vortex street is shown.
99 5.2 Unsteady Cavitating Flows in Convergent-Divergent Nozzles Stutz and Reboud (1997a, 1997b, 2000) have performed a series of experiments to study cavitating flow structure in convergent-divergent nozzles. Time-averaged velocity and vapor volume fraction profiles within the cavity and qualitative description of the cavity behavior have been presented in these studies. Stutz and Reboud (1997b, 2000) have studied unsteady cavitation formed in a different convergent-divergent nozzle. In this particular geometry, the convergent angle is 18 and the divergent part angle is 8. The throat height is 34.3 mm. Cavitation formed in this nozzle is described as unsteady and vapor cloud shedding is typical. The computational domain and the imposed boundary conditions are shown in Figure 5-4. This particular geometry is referred to as Nozzle-1. NoslipNoslipInletOutlet NoslipNoslipInletOutlet Nozzle-2 NoslipNoslipInletOutlet NoslipNoslipInletOutlet Nozzle-1 NoslipNoslipInletOutlet NoslipNoslipInletOutlet Nozzle-2 NoslipNoslipInletOutlet NoslipNoslipInletOutlet Nozzle-1 Figure 5-4. Nozzle-1, computational domain and the imposed boundary conditions. Stutz and Reboud (1997a) have also studied stable sheet cavitation formed in a convergent-divergent nozzle. The convergent part angle is 4.3 and the divergent part angle is 4. The throat height is 43.7 mm. This particular geometry is shown in Figure 5-5 and will be referred to as Nozzle-2.
100 NoslipNoslipInletOutlet Nozzle-1 NoslipNoslipInletOutlet NoslipNoslipInletOutlet Nozzle-2 NoslipNoslipInletOutlet Nozzle-1 NoslipNoslipInletOutlet NoslipNoslipInletOutlet Nozzle-2 Figure 5-5. Nozzle-2, computational domain and the imposed boundary conditions. Reboud et al. (1998) have also conducted a computational study to compare with their experimental findings. In their study, an arbitrary barotropic equation (Delannoy and Kueny 1990, Reboud and Delannoy 1994) is adopted for cavitation modeling. They have indicated that the use of the original kturbulence model leads to a steady-state cavity with no re-entrant jet formed in the closure region, apparently because of high turbulent viscosity induced by the turbulence model. They have modeled the unsteadiness by empirical reduction of the turbulence dissipative terms in the cavitating regions. In the coming sections, cavitating flows forming in the above mentioned nozzles are studied. The experiments indicate unsteady cavitation, shedding clouds of vapor, in Nozzle-1 (Stutz and Reboud 1997b, 2000), and steady sheet cavitation in Nozzle-2 (Stutz and Reboud 1997a). Compared to Nozzle-1, Nozzle-2 has less confinement in the throat section and modest convergent-divergent angles. The original kturbulence model (Jones and Launder 1972) is adopted with no modification. It should be mentioned that the experiments, detailed in Stutz and Reboud (1997a, 1997b, 2000), present the time-averaged velocity and void fraction profiles within the cavity. The unsteadiness observed by them is mostly described qualitatively.
101 5.2.1 Implications of the Speed of Sound Definition In Section 3.4.2, it has been discussed that the speed of sound definition in cavitating regions is important. Basically, any disturbance in the cavitating region propagates at the speed of sound corresponding to the local two-phase mixture. In the incompressible liquid phase, the speed of information propagation is infinite and any disturbance, generated in the cavitating region, is conveyed to the rest of the flow domain. However, the definition of speed of sound in cavitating flows is an open question. Two different approximations, proposed as SoS-1 and SoS-2 in Section 3.4.2, is investigated to shed light on this issue. In particular, SoS-1 is adopted for steady-state computations, because of its good computational stability characteristic. In this section, results obtained from adopting these two definitions are presented and their impact on the time-dependent behavior is discussed. Figure 5-6 shows time evolution of pressure at a point close to the inlet boundary of Nozzle-1. The new interfacial dynamics cavitation model, developed in this study, is used. On the left of Figure 5-6 are the results based on SoS-1, and on the right are the results based on SoS-2. The significant impact of the choice of speed of sound definition is clear. Apparently, SoS-1 does not produce the correct behavior, because the periodic behavior weakens as the cavitation number is decreased. On the other hand, SoS-2, which approximates the fundamental definition of speed of sound, exhibits qualitatively correct behaviors. At higher cavitation numbers, the flow goes to a steady solution after an initial transient, and as the cavitation number is lowered a periodic behavior emerges. The origin of this unsteady behavior will be discussed in more details in the coming section.
102 0 50 100 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 t/T 0 50 100 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 t/T CPCP=1.98=2.25 =1.98=2.25 PPS )1(LSCP =1.98 0 50 100 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 t/T 0 50 100 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 t/T CPCP=1.98=2.25 =1.98=2.25 PPS )1(LSCP =1.98SoS-1SoS-2 0 50 100 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 t/T 0 50 100 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 t/T CPCP=1.98=2.25 =1.98=2.25 PPS )1(LSCP =1.98 0 50 100 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 t/T 0 50 100 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 t/T CPCP=1.98=2.25 =1.98=2.25 PPS )1(LSCP =1.98SoS-1SoS-2 Figure 5-6. The impact of the speed of sound definition on time-dependent behavior. The new interfacial dynamics cavitation model is used. Geometry is Nozzle-1 As discussed in Section 4.2, the steady-state pressure distributions predicted by different cavitation models have been found comparable previously. Differences are mainly in density profiles at the closure region, likely caused by the variations of compressibility characteristics of different cavitation models. The present unsteady simulations support strongly this view as shown in Figure 5-7. This figure compares time-dependent behaviors of the two cavitation models based on the same speed of sound definition, SoS-2. The oscillations of the upstream pressure are considerably different. Model-2 does not produce the unsteady behavior observed in the experiments of Stutz and Reboud (1997b, 2000). The cavity disappears and reappears during its cycle, in
103 contrast to a self oscillating mean cavity observed in the experiments. The results of the new interfacial dynamics cavitation model indicate a stable behavior and as the cavitation number is decreased a periodic behavior with mean cavity emerges. PPS PPS 0 50 100 -1 -0.5 0 0.5 1 1.5 0 50 100 -1 -0.5 0 0.5 1 1.5 CPCPt/Tt/T=2.25 =2.25 PPS PPSSoS-2SoS-2 0 PPS PPS 0 50 100 -1 -0.5 0 0.5 1 1.5 0 50 100 -1 -0.5 0 0.5 1 1.5 CPCPt/Tt/T=2.25 =2.25 PPS PPSSoS-2SoS-2 0 Figure 5-7. Comparison of the time-dependent behavior of different cavitation models. Geometry is Nozzle-2. The speed of sound definition has been investigated for computations of unsteady cavitation formed in Nozzle-1, which has been observed in the experiments of Stutz and Reboud (1997b, 2000). As already mentioned, Stutz and Reboud (1997a) have reported steady-state cavity dynamics in their experiments on sheet cavitation formed in Nozzle-2. shows the time evolution of pressure at a point close to the inlet boundary of Nozzle-2. Both the new interfacial dynamics cavitation model and Model-2 are used. As Figure 5-8
104 seen on the right plot of , both models exhibit unsteady behavior if SoS-1 is adopted. This is obviously not consistent with the essentially steady-state cavity observed in the experiments (Stutz and Reboud 1997a). But if SoS-2 is adopted, which is an approximation to the fundamental speed of sound, the correct stable behavior is captured by both cavitation models. Figure 5-8 Figure 5-8. The impact of the speed of sound definition on time-dependent behavior. Geometry is Nozzle-2. The cavitation number is 0.67. -0.65 -0.6 -0.55 -0.5 -0.45 -0.4 t/TCp 0 50 100 -0.7 -0.65 -0.6 -0.55 -0.5 -0.45 -0.4 t/TCp New Model Model-2 Model-2 )1(LSCP PPSSoS-2SoS-1 -0.65 -0.6 -0.55 -0.5 -0.45 -0.4 t/TCp 0 50 100 -0.7 -0.65 -0.6 -0.55 -0.5 -0.45 -0.4 t/TCp Model-2 Model-2 New Model )1(LSCP )1(LSCP PPS PPSSoS-2SoS-1 -0.65 -0.6 -0.55 -0.5 -0.45 -0.4 t/TCp 0 50 100 -0.7 -0.65 -0.6 -0.55 -0.5 -0.45 -0.4 t/TCp New Model Model-2 Model-2 )1(LSCP PPSSoS-2SoS-1 -0.65 -0.6 -0.55 -0.5 -0.45 -0.4 t/TCp 0 50 100 -0.7 -0.65 -0.6 -0.55 -0.5 -0.45 -0.4 t/TCp Model-2 Model-2 New Model )1(LSCP )1(LSCP PPS PPSSoS-2SoS-1 Based on these numerical tests, it seems that SoS-2, gives qualitatively correct time-dependent behaviors. This is not surprising because SoS-2 is an approximation to the fundamental definition of speed of sound. The mean flow direction is chosen as an approximation for the path to compute the partial derivative in the fundamental definition. Such an approximation is needed because the actual path is not known.
105 This test study, regarding the speed of sound definition and cavitation models, has shown that the new interfacial dynamics cavitation model along with SoS-2 definition is more likely to produce the correct unsteady behavior and is further employed to investigate the unsteady cavitation. 5.2.2 Cavitation Instability (Auto-Oscillations) The results presented in the previous section have shown that a periodic behavior emerges in Nozzle-1 as the cavitation number is lowered, whereas the cavitation appears to be stable in Nozzle-2. In this section the goal is to identify the origin of the cavitation instability observed in Nozzle-1 and to characterize it. Callenaera et al. (2001) have proposed two types of instability for an attached cavity: intrinsic instabilities and system instabilities. An intrinsic instability originates within the cavity itself, such as the re-entrant jet instability causing the shedding of cloud cavitation. On the other hand, system instability is because of the interaction between the cavity and other components of the hydraulic system such as the inlet and outlet lines. Watanabe et al. (1998) and Franc (2001) have discussed that cavitation surge is of system instability type, for which the dynamic behavior of the cavity is strongly coupled to its environment. With a one dimensional model it has been shown that fluctuations in the incoming mass flow rate causes pressure oscillations in the hydraulic system. Duttweiler and Brennen (2002) have detailed experiments of cavitation surge instability on a cavitating propeller in a water tunnel. It is reported that the onset of instability is observed at decreasing cavitation numbers. The periodic pressure oscillations resulting from this instability have been measured away from the cavitating region.
106 In light of these studies, the findings obtained for the Nozzle-1 case are discussed. shows the shapes of the cavity during the auto-oscillation cycle. The computations have captured a mean cavity and recirculation zone. Snapshot of the unsteady cavitation (Stutz and Reboud 1997b) is also included to assess the present findings. The main cavity body is in good agreement with the experiment, as shown in . The experimental photograph indicates a large scale structure that rolls up and sheds clouds of vapor, which is not captured in the present computations. Kubota et al. (1989) have associated the cloud cavitation motion with the large scale velocity. The reason that the cloud cavitation is not captured in the present results may be due to the adoption of the original kturbulence model (Jones and Launder 1974), which induces high turbulent viscosities in regions of turbulence production. The typical Strouhal number of the shedding of cloud cavitation is found to be 0.30 in various studies (Stutz and Reboud 1997b, Callenaere et al. 2001). In the present results, no shedding of cloud cavitation is captured; hence, the Strouhal number is not in agreement with the above value. The reduced frequency of the auto-oscillations of the cavity, fc/U, based on characteristic length and the reference velocity is computed as 0.08 based on the cavity cycle shown in Figure 5-9. This number is in agreement with the reduced frequency (0.07) of the cavitation surge instability on a propeller, presented in Duttweiler and Brennen (2002), where no periodic shedding of cloud cavitation is observed. It appears that there are two dominant scales related to the observed unsteadiness, large scale motion responsible for vortex shedding and cloud cavitation, and smaller scale motion for attached cavity. The Strouhal numbers associated with these two regimes are different in Figure 5-9Figure 5-9
107 magnitudes. The present modeling strategy seems to capture the small scale transient behavior more easily. Figure 5-10 shows the periodic oscillation of the upstream pressure in Nozzle-1. The similar behavior of the far field pressure is reported as cavitation surge in Duttweiler and Brennen (2002) for a cavitating propeller in a water tunnel. As previously shown in the right plot of , this periodic behavior emerges as the cavitation number is decreased. Figure 5-6 Callenaere et al. (2001) have studied the effect of confinement in a venturi type test section and have found that the auto oscillations develop as the confinement is increased, indicating the prominent role of pressure gradient along the channel. In their study, the pressure gradient has also been controlled by varying the divergent angle of the test section. The auto-oscillation disappears as the divergent angle is decreased. If the results obtained from two nozzles (Nozzle-1 and Nozzle-2) are investigated in a collective approach, the prominent role of the pressure gradient in generating the auto-oscillations is clearly identified. To summarize qualitatively, Nozzle-1 has a narrow throat section with large convergent-divergent angles, whereas Nozzle-2 has a wider throat section with smaller convergent-divergent angles; hence, the flow is less confined in the case of Nozzle-2. In Nozzle-1, the cavity incepts and grows in length and thickness contracting the effective area at the throat section, which leads to auto oscillation of the cavity. However in Nozzle-2 steady-state sheet cavitation is observed. Such characteristics of the present numerical simulations are consistent with the experimental observations of Callenaere et al (2001).
108 Snapshot of unsteady cavitation Stutz and Reboud (1997b)Present computationt/T=0.0t/T=7.0t/T=10.0t/T=12.6 Snapshot of unsteady cavitation Stutz and Reboud (1997b)Present computationt/T=0.0t/T=7.0t/T=10.0t/T=12.6 Figure 5-9. Cavity shapes and the recirculation zone during the instability cycle. Geometry is Nozzle-1. The new interfacial dynamics cavitation model is used. The cavitation number is 1.98 based on time-averaged quantities. The experimental photo is adopted from Stutz and Reboud (1997b) with the permission of Springer-Verlag.
109 0 20 40 60 80 100 120 140 160 180 200 -0.5 -0.45 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 t/TCp Figure 5-10. The periodic oscillation of the upstream pressure in Nozzle-2. The cavitation number is 1.98 based on time-averaged quantities. 5.2.3 Two-Phase Flow Structure Time-averaged velocity and vapor volume fraction profiles are compared with experimental data of Stutz and Reboud (1997a, 2000) to assess the present computational results of Nozzle-1 and Nozzle-2 cases. Figure 5-11 and show the time-averaged velocity and vapor volume fraction profiles within cavity formed in Nozzle-1. The boundary of the cavitating region is also included in Figure 5-11. The computations capture the main cavity body. The total time-averaged cavity boundary is not captured because the shedding of cloud cavitation does not occur in the computations due to the high turbulent viscosity induced by the original kturbulence model. Reasonable agreement is observed in the velocity profiles Figure 5-12
110 given in Figure 5-11, especially in the core of the reverse flow. The vapor volume fraction profiles obtained in the experiments shows high mixture content, which is typical in cloud cavitation. As expected, the computations do not match the experimental data quantitatively. Nevertheless the overall trends are similar. 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 14 16 18 20 X(mm)Y(mm)Time-averaged velocity profiles within the cavity 10.00.00.00.00.010.010.010.0 Present computation Exp. data Boundary of the Cavitating flow (Exp.) 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 14 16 18 20 X(mm)Y(mm)Time-averaged velocity profiles within the cavity 10.00.00.00.00.010.010.010.0 Present computation Exp. data Boundary of the Cavitating flow (Exp.) Figure 5-11. Time averaged velocity profiles (u/U ref ) within the cavity formed in Nozzle-1. The vertical scale is the distance from the wall. Experimental data is from Stutz and Reboud (2000). The cavitation number is 1.98. Figure 5-13 and shows the velocity and vapor volume fraction profiles, respectively, within the cavity formed in Nozzle-1. Recall that a steady-state computation has been performed for this particular geometry and the corresponding results are documented in Figure 4-30 and Figure 4-31. By comparison, the time-averaged results have improved the agreement with the time-averaged experimental data. Considering that Figure 5-14
111 all the experimental data are measurements within the cavity, the computations give a good agreement in most part of the cavity. The agreement in the reverse flow is improved as compared to the results based on steady-state formulation. 10 20 30 40 50 60 70 80 0 2 4 6 8 10 12 14 16 18 Time-averaged vapor-volume fraction profilesX(mm)Y(mm) 0.0252525250.00.00.0V(%) Present computationExp. data 10 20 30 40 50 60 70 80 0 2 4 6 8 10 12 14 16 18 Time-averaged vapor-volume fraction profilesX(mm)Y(mm) 0.0252525250.00.00.0V(%) 10 20 30 40 50 60 70 80 0 2 4 6 8 10 12 14 16 18 Time-averaged vapor-volume fraction profilesX(mm)Y(mm) 0.0252525250.00.00.0V(%) Present computationExp. data Figure 5-12. Vapor volume fraction profiles within the cavity formed in Nozzle-1. The vertical scale is the distance from the wall. Experimental data is from Stutz and Reboud (2000).The cavitation number is 1.98. The vapor volume fraction profiles have also improved by adopting unsteady computations. The disagreement is in the reverse flow region. The computations indicate mostly liquid content, whereas the experimental results show a mixture flow. This also explains, in part, the disagreement of the velocity profiles in the reverse flow. The reverse mixture flow, observed in the experiments, has higher velocity due to its lower density. However, the high density slows down the motion of the reverse flow, captured in the computations.
112 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 0.12 X/Lcav(%)Y/LcavTime-averaged velocity profiles within the cavity 010 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 0.12 X/Lcav(%)Y/LcavTime-averaged velocity profiles within the cavity 010 Figure 5-13. Time-averaged velocity profiles (u/U ref ) within the cavity formed in Nozzle-2. The vertical and horizontal scales are normalized by Lcav=80mm to be consistent with the experimental data. The vertical scale is the distance from the wall. Experimental data is from Stutz and Reboud (1997a). The cavitation number is 0.67.
113 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 X/Lcav(%)Y/LcavTime-averaged vapor volume fraction profiles V (%)0100 Exp. dataPresent computation 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 X/Lcav(%)Y/LcavTime-averaged vapor volume fraction profiles V (%)0100 Exp. dataPresent computation Figure 5-14. Time-averaged vapor volume fraction profiles within the cavity formed in Nozzle-2. The vertical and horizontal scales are normalized by L cav =80mm to be consistent with the experimental data. The vertical scale is the distance from the wall. Experimental data is from Stutz and Reboud (1997a). The cavitation number is 0.67.
CHAPTER 6 CONCLUSIONS AND FUTURE DIRECTIONS 6.1 Conclusions A comprehensive computational methodology is developed for turbulent cavitating flows. The computational capability is applied to axisymmetric projectiles and a planar hydrofoil for steady-state simulation and to convergent-divergent nozzles for both steady-state and time-dependent simulation. The ensemble-averaged Navier-Stokes equations, cast in their conservation form, along with a volume fraction transport equation, are employed. The flow is computed in both phases with the vapor pressure recovered inside the cavity via mass transfer. The two-equation kmodel is adopted for turbulence closure. To ensure convergent and stable computations of turbulent cavitating flows, a unified incompressible-compressible pressure-based algorithm is established for steady-state simulations as the first step. A pressure-density coupling scheme is developed to handle the large density jump, O(1000), at the phase boundary. The resulting pressure-correction equation shares common features with that of high-speed flows, exhibiting a convective-diffusive type, instead of only diffusive type. Specifically, the convective effect in the pressure-correction equation is important as the vapor volume fraction increases in the cavitating region; elsewhere the pressure-correction equation returns to a purely diffusive type. Similar to high-speed cases, upwinded density interpolation in mass flux 114 computations also aids convergence of the cavitating flow computations. Several issues
115 related to grid resolution, convection schemes, and turbulence modeling are investigated. The recirculation zone in the closure region is better resolved with grid refinement. While satisfactory agreement between numerical solution and experimental measurements in wall pressure distribution can be maintained with variations in grid resolution and parameters in the empirical cavitation model, other aspects such as the density distribution are found to exhibit higher sensitivity to them. A remark has been made regarding the potential of transport equation-based cavitation models to capture the baroclinic vorticity generation in the closure region of cavity. An analysis of the mass and momentum conservation at a phase-change liquid-vapor interface is conducted in the context of homogeneous equilibrium flow theory. This analysis has shed light into the source of empiricism in existing transport equation-based cavitation models. Based on this analysis of interfacial dynamics, a transport equation-based interfacial dynamics cavitation model is developed. The new model replaces the empiricism by computing the interfacial velocities and it is general for time-dependent computations. The unified pressure-based algorithm is further employed to assess three versions of the transport equation-based empirical cavitation models in literature and the new interfacial dynamics cavitation model. The steady-state computations have shown that comparable pressure distributions are produced by all the cavitation models considered. Noticeable variations have been observed in density distribution, implying that each cavitation model has different compressibility characteristics. The new interfacial dynamics cavitation model is further applied to convergent-divergent nozzles to study the two-phase flow structure of sheet cavitation. It is found that the baroclinic vorticity
116 generation is important for the turbulence production in the cavity closure region. The abrupt variation of density has certain implications on the near-wall turbulence treatment with the kturbulence model. Due to depression of density in the cavitating regions, the local Reynolds number reduces substantially and the essence of the equilibrium turbulence assumption of the wall function formulation may be lost in certain cases. The issue weakens the convergence and stability of cavitating flow computations. The problem is identified and overcome by adopting wall-function consistent grids. For time-dependent simulations, the ideas embedded in the unified pressure-based algorithm are further employed in the pressure-based operator splitting algorithm. The sequence of the basic operations in the pressure-based operator splitting algorithm are detailed for turbulent cavitating flow computations. The algorithm is enhanced for time-accurate computations by repeating the corrector steps. Proper coupling of the cavitation and turbulence models into the corrector steps is important as it affects the stability of the computation. The implications of the local speed of sound definition in the two-phase mixture are investigated for time-dependent computations. Two approximations regarding the speed of sound definition in the two-phase mixture have been proposed. It is shown that the definition of local speed of sound of the two-phase mixture has significant impact on the unsteady behavior captured in the simulations. The formula that approximates the fundamental definition of speed of sound is found to produce the experimentally observed unsteady behavior. Cavitation instability is observed as auto-oscillations of the attached cavity formed in the convergent-divergent nozzle. The upstream pressure oscillates periodically in response to cavity oscillations. The reduced frequency of the auto oscillations is found to
117 be in good agreement with the experimental value. The effect of pressure gradient on the cavitation instability is highlighted by considering two convergent-divergent nozzles with different degree of geometric confinement around the cavitation region. To conclude, the computational methodology has exhibited capabilities to predict both qualitative and quantitative features of turbulent cavitating flows. The simulations have helped get reasonable insight into the complex flow structure of cavitation. 6.2 Future Directions The present research can be extended in the following directions: Turbulence modeling for unsteady two-phase mixtures. Development of interface tracking capability that will enable the information of interface velocity to be used in the new interfacial dynamics cavitation model. Coupling of the energy equation to consider compressible flow situations and as well as to study cavitation of cryogenic liquids. Application to rotational flows encountered in turbomachinery analysis. Further investigation of the computation of local speed of sound in two-phase mixtures.
118 Collier, J.G., and Thome, J.R., 1994, Convective Boiling and Condensation, 3 LIST OF REFERENCES Ahuja, V., Hosangadi, A., and Arunajatesan, S., 2001, â€œSimulations of Cavitating Flows Using Hybrid Unstructured Meshes,â€ J. Fluid Eng-T. ASME, 123, pp.331-340. Arakeri, V.H., 1975, â€œViscous Effects on the Position of Cavitation Separation from Smooth Bodies,â€ Arakeri, V.H., and Acosta, A.J., 1973, â€œViscous Effects in the Inception of Cavitation on Axisymmetric Bodies,â€ J. Fluid Eng-T. ASME, 95(4), pp. 519-527. Batchelor, G.K., 1967, An Introduction to Fluid Dynamics, Cambridge University Press, New York. Brackbill, J.U., Kothe, D.B. and Zemach, C., 1992, â€œA Continuum Method for Modeling Surface Tension,â€ J. Comput. Phys., 100, pp. 335-354. Brennen, C.E., 1995, Cavitation and Bubble Dynamics, Oxford University Press, New York. Bressloff, N.W., 2001, â€œA Parallel Pressure Implicit Splitting of Operators Algorithm Applied to Flows at All Speeds,â€ Int. J. Numer. Meth. Fl., 36, pp. 497-518. Callenaere, M., Franc J.P., Michel, J.M., and Riondet M., 2001, â€œThe Cavitation Instability Induced by the Development of a Re-entrant Jet,â€ J. Fluid Mech., 444, pp.223-256. Carey, V.P., 1992, Liquid-Vapor Phase-Change Phenomena, Hemisphere, Washington DC. Ceccio, S., Gowing, S., and Shen, Y.T., 1997, â€œThe Effects of Salt Water on Bubble Cavitation,â€ J. Fluid Eng-T. ASME, 119, pp. 155-163. Chen, Y., and Heister, S.D., 1994, â€œA Numerical Treatment for Attached Cavitation,â€ J. Fluid Eng-T. ASME, 116, pp. 613-618. Chen, Y., and Heister, S.D., 1996, â€œModeling Hydrodynamics Nonequilibrium in Cavitating Flows,â€ J. Fluid Eng-T. ASME, 118, pp. 172-178. n, Carendon Press, Oxford. rd editio
119 Delannoy, Y., and Kueny, J.L., 1990, â€œCavity Flow Predictions Based on the Euler Equations,â€ ASME Cavitation and Multiphase Flow Forum. Deshpande, M., Feng, J., and Merkle, C.L., 1997, â€œNumerical Modeling of the Thermodynamic Effects of Cavitation,â€ J. Fluid Eng-T. ASME, 119, pp. 420-427. Duttweiler, M.E., and Brennen, C.E., 2002, â€œSurge Instability on a Cavitating Propeller,â€ J. Fluid Mech. 458, pp. 133-152. Edwards, J.R., Franklin, R.K., and Liou, M.S., 2000, â€œLow-Diffusion Flux-Splitting Methods for Real Fluid Flows with Phase Transitions,â€ AIAA J., 38(9), pp. 1624-1633. Ferziger, J.H., and Peric M., 1999, Computational Methods for Fluid Dynamics, 2 nd edition, Springer Verlag, Berlin, Germany. Franc, J.P., 2001, â€œPartial Cavity Instabilities and Re-Entrant Jet,â€ CAV2001, Proc. 4th International Symposium on Cavitation, CAV2001:lecture.002, California Institute of Technology, Pasadena, CA, ( http://cav2001.library.caltech.edu ). Franc, J.P., Avellan, F., Belahadji, B., Billard, J.Y., Briancon, L., Marjollet, Frechou, D., Fruman, D. H., Karimi, A., Kueny, J.L., and Michel, J. M., 1995, La Cavitation: Mecanismes Physiques et Aspects Industriels, Presses Universitaires de Grenoble, Grenoble (EDP Sciences), pp. 28, 141, 142, 144, 449, France. Franc, J.P., and Michel, J.M., 1985, â€œAttached Cavitation and the Boundary Layer: Experimental Investigation and Numerical Treatment,â€ 154, pp.63-90. Francois, M., 1998, â€œA Study of the Volume of Fluid Method for Moving Boundary Problems,â€ M.Sc. Thesis, Embry Riddle Aeronautical University, Daytona Beach FL, USA. Gopalan, S., and Katz, J., 2000, â€œFlow Structure and Modeling Issues in the Closure Region of Attached Cavitation,â€ Phys. Fluids, 12(4), pp. 895-911. Harten, A., 1984, â€œHigh Resolution Schemes for Hyberbolic Conservation Laws,â€ J. Comput. Phys., 49, pp. 357-393. He, X., Senocak, I., Shyy, W., Gangadharan, S.N., and Thakur, S., 2000, â€œEvaluation of Laminar-Turbulent Transition and Equilibrium Near-wall Turbulence Models,â€ Numer. Heat Tr. A-Appl, 37, pp. 101-112. Incropera, F.P., and Dewitt, D.P., 1996, Fundamentals of Heat and Mass Transfer, 4 th edition, John Wiley & Sons, New York.
120 Issa, R.I., 1985, â€œSolution of the Implicitly Discretized Fluid Flow Equations by Operator-Splitting,â€ J. Comput. Phys., 62, pp. 40-65. Johnson, R.W., 1998, The Handbook of Fluid Dynamics, CRC Press, Boca Raton, FL. Jones, W.P., and Launder, B.E., 1972, â€œThe Prediction of Laminarization with a Two-Equation Model of Turbulence,â€ Int. J. Heat Mass Tran., 15, pp. 301-314. Joseph, D.D., 1995 â€œCavitation in a Flowing Liquid,â€ Phys. Review E, 51(3), pp. 1649. Karki, K.C., and Patankar S.V., 1989, â€œPressure Based Calculation Procedure for Viscous Flows at All Speeds in Arbitrary Configurations,â€ AIAA J., 27(9), pp. 1167-1174. Kato, H., Yamaguchi, H., and Maeda, M., 1988, â€œDirect Measurements of Shearing Stress an Heat Transfer on a Flat Plate Covered with a Sheet Cavity,â€ Proceedings of ASME Cavitation and Multiphase Flow Forum, Cincinnati, OH. Katz, J., 1984, â€œCavitation Phenomena within Regions of Flow Separation,â€ J. Fluid Mech., 140, pp. 397-436. Kawanami, Y., Kato H., Tanimura, M., and Tagaya, Y., 1997, â€œMechanism and Control of Cloud Cavitation,â€ J. Fluid Eng-T. ASME, 119, pp. 788-794. Kinnas, S.A., and Mazel, C.H., 1993, â€œNumerical versus Experimental Cavitation Tunnel,â€ J. Fluid Eng-T. ASME, 115, pp. 760-765. Kirschner, I.N., 2001, â€œResults of Selected Experiments Involving Supercavitating Flows,â€ VKI Lecture Series on Supercavitating Flows, VKI Press, Brussels, Belgium. Knapp, R. T., Daily, J.W. and Hammitt, F. G., 1970, Cavitation, McGraw-Hill, New York. Kubota, A., Kato, H., and Yamaguchi, H., 1992, â€œA New Modelling of Cavitating Flows: A Numerical Study of Unsteady Cavitation on a Hydrofoil Section,â€ J. Fluid Mech., 240, pp.59-96. Kubota, A., Kato, H., Yamaguchi, H., and Maeda, M., 1989, â€œUnsteady Structure Measurement of Cloud Cavitation on a Foil Section Using Conditional Sampling Technique,â€ J. Fluid Eng-T. ASME, 111, pp. 204-210. Kuhn de Chizelle, Y., Ceccio, S.L., and Brennen, C.E., 1995, â€œObservations and Scaling of Traveling Bubble Cavitation,â€ J. Fluid Mech., 293, pp.96-126. Kunz, R.F., Boger, D. A., Stinebring, D. R., Chyczewski, T.S., Lindau, J.W., and Gibeling, H. J., 2000 â€œA Preconditioned Navier-Stokes Method for Two-phase Flows with Application to Cavitation,â€ Comput Fluids, 29, pp.849-875.
121 Launder, B.E., and Spalding, D. B., 1974, â€œThe Numerical Computation of Turbulent Flows,â€ Comput. Method. Appl. M., 3, pp. 269-289. Lecoffre, Y., 1999, Cavitation Bubble Trackers, A.A. Balkema, Brookfield, VT. Lesieur, M., 1993, Turbulence in Fluids, Second revised edition, Kluwer Academic Publishers, Dordrecht, The Netherlands. Li, C.Y., Ceccio, S.L., 1996, â€œInteraction of Single Travelling Bubbles with the Boundary Layer and Attached Cavitation,â€ J. Fluid Mech., 322, pp. 329-353. Merkle, C.L., Feng, J., and Buelow, P.E.O., 1998 â€œComputational Modeling of the Dynamics of Sheet Cavitation,â€ Proc. 3 rd International Symposium on Cavitation, Grenoble, France. Moukalled, F., and Darwish, M., 2000, â€œA Unified Formulation of the Segregated Class of Algorithms for Fluid Flow at All Speeds,â€ Numer. Heat Tr. B-Fund, 37, pp.103-139. Oliveira, P.J., and Issa, R.I., 2001, â€œAn Improved PISO Algorithm for the Computation of Buoyancy-Driven Flows,â€ Numer. Heat Tr. B-Fund, 40, pp.473-493. Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington DC. Patankar, S.V., and Spalding, D.B., 1972, â€œA Calculation Procedure for Heat, Mass and Momentum Transfer in Three Dimensional Parabolic Flows,â€ Int. J. Heat Mass Tran, 15, pp. 1787-1806. Pope, S.B., 2000, Turbulent Flows, Cambridge University Press, Cambridge, UK. Rayleigh L., 1917, â€œOn the Pressure Developed in a Liquid During the Collapse of a Spherical Cavity,â€ Phil. Mag., 34, pp. 94-98. Reboud, J.L., and Delannoy, Y., 1994, â€œTwo-Phase Flow Modelling of Unsteady Cavitation,â€ Proc. 2 nd International Symposium on Cavitation, Tokyo, Japan. Reboud, J.L., Stutz, B., and Coutier, O., 1998, â€œTwo-Phase Flow Structure of Cavitation: Experiment and Modelling of Unsteady Effects,â€ Proc. 3 rd International Symposium on Cavitation, Grenoble, France. Reisman, G.E., Wang, Y.C., and Brennen, C.E., 1998, â€œObservations of Shock Waves in Cloud Cavitation,â€ J. Fluid Mech., 355, pp. 255-283.
122 Reynolds, W.C., 1979, Thermodynamic Properties in SI, Stanford University, Palo Alto, CA. Rhie, C., and Chow, W., 1983, â€œNumerical Study of the Turbulent Flow Past an Airfoil with Trailing Edge Separation,â€ AIAA J., 21, pp. 1525-1532. Rouse, H. and McNown, J.S., 1948, â€œCavitation and Pressure Distribution, Head Forms at Zero Angle of Yaw,â€ Studies in Engineering, Bulletin 32, State University of Iowa Aemes, Iowa. Schetz, J.A., 1984, Foundations of Boundary Layer Theory for Momentum, Heat, and Mass Transfer, Prentice-Hall, Englewood Cliffs, NJ. Schlichting, H., 1979, Boundary Layer Theory, McGraw-Hill, New York. Senocak, I., and Shyy, W., 2001, â€œNumerical Simulation of Turbulent Flows with Sheet Cavitation,â€ CAV2001, Proc. 4th International Symposium on Cavitation, CAV2001A7.002, California Institute of Technology, Pasadena, CA, ( http://cav2001.library.caltech.edu ). Senocak, I., and Shyy, W., 2002, â€œA Pressure-Based Method for Turbulent Cavitating Flows,â€ J. Comput. Phys., 176, pp. 363-383. Shen, Y., and Dimotakis, P., 1989 â€œThe Influence of Surface Cavitation on Hydrodynamic Forces,â€ Proc. 22 nd ATTC, St. Johns, pp. 44-53. Shyy, W., and Braaten, M.E., 1988, â€œAdaptive Grid Computation for Inviscid Compressible Flows Using a Pressure Correction Method,â€ Proc. AIAA/ASME/SIAM/APS First National Fluid Dynamics Congress, AIAA-CP 888, pp. 112-120. Shyy, W., 1994, Computational Modeling for Fluid Flow and Interfacial Transport, Elsevier, Amsterdam, The Netherlands, (revised printing 1997). Shyy, W., and Thakur, S.S., 1994a, â€œA Controlled Variation Scheme in a Sequential Solver for Recirculating Flows. Part I. Theory and Formulation,â€ Numer. Heat Tr.B-Fund., 25(3), pp. 245-272. Shyy, W., and Thakur, S.S., 1994b, â€œA Controlled Variation Scheme in a Sequential Solver for Recirculating Flows. Part II. Applications,â€ Numer. Heat Tr.A-Appl., 25(3), pp. 273-286. Shyy, W., Thakur, S.S., Ouyang, H., Liu, J. and Blosch, E., 1997, Computational Techniques for Complex Transport Phenomena, Cambridge University Press, New York.
123 Singhal, A.K., Li, H., Athavale, M.M., and Jiang, Y., 2001, â€œMathematical Basis and Validation of the Full Cavitation Model,â€ ASME Paper FEDSM2001-18015, Proc. of 2001 ASME Fluids Engineering Division Summer Meeting. Singhal, A.K., Vaidya, N., and Leonard, A.D., 1997, â€œMulti-dimensional Simulation of Cavitating Flows Using a PDF Model for Phase Change,â€ ASME Paper FEDSM97-3272, Proc. of 1997 ASME Fluids Engineering Division Summer Meeting. Spalart, P.R., 2000, â€œStrategies for Turbulence Modelling and Simulations,â€ Int. J. Heat Fluid Fl., 21, pp.252-263. Stinebring, D.R., Billet, M.L., Lindau, J.W., and Kunz, R.F., 2001, â€œDeveloped Cavitation-Cavity Dynamics,â€ VKI Lecture Series on Supercavitating Flows, VKI Press, Brussels, Belgium. Stutz, B., and Reboud, J.L., 1997a, â€œTwo-phase Flow Structure of Sheet Cavitation,â€ Phys. Fluids, 9(12), pp. 3678-3686. Stutz, B., and Reboud, J.L., 1997b, â€œExperiments on Unsteady Cavitation,â€ Exp. Fluids, 22, pp. 191-198. Stutz, B., and Reboud, J.L., 2000, â€œMeasurements within Unsteady Cavitation,â€ Exp. Fluids, 29, pp. 545-552. Tannehill, J.C., Anderson, D.A., Pletcher, R.H., 1997, Computational Fluid Mechanics and Heat Transfer, Second edition, Taylor & Francis, Washington, DC, USA. Thakur, S.S., and Wright, J.F., 2002, â€œAn Operator-Splitting Algorithm for Unsteady Flows at All Speeds in Complex Geometries,â€ submitted for publication. Thakur, S.S., Wright, J.F., Shyy, W., and Udaykumar, H., 1997, â€œSEAL: A Computational Fluid Dynamics and Heat Transfer Code for Complex 3-D Geometries,â€ University of Florida, Gainesville. Vandoormal, J.P., and Raithby, G.D., 1984, â€œEnhancements of the SIMPLE Method for Predicting Incompressible Fluid Flows,â€ Numer. Heat Tr., 7, pp. 147-163. Venkateswaran, S., Lindau, J.W., Kunz, R.F., and Merkle, C.L., 2001, â€œPreconditioning Algorithms for the Computation of Multi-phase Mixtures Flows,â€ AIAA-2001-0125, AIAA 39 th Aerospace Sciences Meeting and Exhibit. Ventikos, Y., and Tzabiras, G., 2000, â€œA numerical Method for the Simulation of Steady and Unsteady Cavitating Flows,â€ Comput Fluids, 29, 63-88. Versteeg, H.K., and Malalasekera W., 1995, An Introduction to Computational Fluid Dynamics, Longman, London, UK.
124 Yoshihara, K., Kato, H., and Yamaguchi, H., 1988, â€œExperimental Study on the Internal Flow of a Sheet Cavity,â€ ASME Cavitation and Multiphase Flow Forum, Cincinnati, OH, pp.114-11. Wallis, G.B., 1969, One-dimensional Two-Phase Flow, McGraw Hill, New York City, NY. Wang, G., Senocak, I., Shyy, W., Ikohagi, T., Cao, S., 2001, â€œDynamics of Attached Turbulent Cavitating Flows,â€ Prog. Aerosp. Sci., 37, pp.551-581. Watanabe, S., Tsujimoto, Y., Franc, J.P., and Michel, J.M., 1998, â€œLinear Analyses of Cavitation Instabilities,â€ Proc. 3 rd International Symposium on Cavitation, Grenoble, France. White, F.M., 1974, Viscous Fluid Flow, McGraw-Hill, New York. Wilcox, D.C., 1993, Turbulence Modeling for CFD, DCW Industries, La Canada, California. Wilcox, D.C., 1997, Basic Fluid Mechanics, DCW Industries, La Canada, California. Williamson, C.H.K., 1989, â€œOblique and Parallel Mode of Vortex Shedding in the Wake of a Cylinder at Low Reynolds Number,â€ J. Fluid Mech. 206, pp.579-627.
BIOGRAPHICAL SKETCH nan enocak was born in Ankara, Turkey, in May 1976. He received his B.S. degree in Mechanical Engineering from the Middle East Technical University, Ankara, Turkey in June 1998. Soon after receiving his degree, he enrolled in the Ph.D. program in Aerospace Engineering at the University of Florida, Gainesville, USA in August 1998. 125