UFDC Home  myUFDC Home  Help 



Full Text  
FILTERBASED MODELING OF UNSTEADY TURBULENT CAVITATING FLOW COMPUTATIONS By JIONGYANG WU A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2005 Copyright 2005 by Jiongyang Wu To my Parents; my wife, Weishu; and my son, Andy, for their love and support ACKNOWLEDGMENTS I would like to express my sincere thanks and appreciation to my advisor Professor Wei Shyy. His thorough guidance and support benefited me tremendously in exploring and pursuing my research interest. I thank him for his boundless patience and generous attitude, and his enduring enthusiasm in educating me both in research and personal development. I have benefited much from the collaboration with Professor Stein T. Johansen. I would like to thank Professor Renwei Mei for providing some helpful suggestions and serving on my committee. Additionally, I would like to thank Professors Louis N. Cattafesta, David W. Mikolaitis, and Jacob N. Chung for serving on my committee, and Dr. Siddharth Thakur for generously sharing his experience. I have had the privilege to work with all the individuals of the Computational ThermoFluids group, to whom I also give my thanks. I would like to deeply thank my parents. Their encouragement, support, trust, and love have given me power and strength through these years. My wife, Weishu Bu, has been with me and supported me all the time. My son, Andy Wu, is another source of my invaluable wealth. No words can possibly express my gratitude and love to them. Finally, I would like to thank the NASA Constellation University Institute Program (CUIP) for financial support. TABLE OF CONTENTS page A C K N O W L E D G M E N T S ................................................................................................. iv LIST OF TABLES ............. ................... ......... .... ... ............... .. vii LIST O F FIG U R E S ......................................................... ......... .. ............. viii L IST O F SY M B O L S .... ....................................................... .. ....... .............. xi ABSTRACT ........ .............. ............. ...... ...................... xiv CHAPTER 1 IN TR OD U CTION ............................................... .. ......................... .. 1.1 B background of C aviation ............................................................................ 1 1.1.1 C aviation Types in Fluids....................................................... ............ .. 1.1.2 Cavitation Inception and Param eter ........................................ ..................2 1.2 Research Motivations and Objectives .......................................... ...............3 1.2.1 Challenges and M otivations ........................................ ...... ............... 3 1.2 .2 R research O bjectiv es ........................................................................ ...... 8 2 NAVIERSTOKES EQUATIONS AND TURBULENCE MODELS ....................9 2 .1 G ov earning E qu action s .................................................................... .............. .9 2.2 Turbulence M odel ....... .... .................................. ........ .. .. .. ............ 11 2.3 FilterBased RAN S M odel .............................................................................12 2 .3.1 L literature R review ............................................................ ....................14 2.3.2 FilterBased M odel Concept (FBM ) .................................. ............... 16 2.3.3 Modeling Implementation ............................ ..............17 2.4 Assessing FBM in TimeDependent SinglePhase Flow............... ..................18 3 CA V ITA TION M OD ELS ................................................ .............................. 32 3.1 Governing Equations ....................... ............................ 32 3.2 Literature Review of Cavitation Studies.................................... ............... 33 3.2.1 Cavitation Com pressibility Studies ......................... ....................... .... 35 3.2.2 Cavitation Studies on Fluid Machinery Components and Systems............38 3.3 Transport Equationbased Cavitation Model (TEM)................ .............. ....41 4 NUMERICAL METHODS .................................... ............................. 45 4.1 PressureB ased A lgorithm ............................................................. ................ .... 46 4.2 Pressure Implicit Splitting of Operators (PISO) Algorithm for Unsteady Com putations ............................ ..................... 47 4.3 SpeedofSound (SoS) Numerical Modeling............................................50 5 ASSESSING TIMEDEPENDENT TURBULENT CAVITATION MODELS .......53 5.1 Cavitating Flow through a HollowJet Valve ............... ............. ....................53 5.1.1 Steady and Unsteady Turbulent Cavitating Flows ...................................55 5.1.2 Impact of Speedofsound Modeling..................................... .................60 5.2 Turbulent Cavitating Flow through a ConvergentDivergent Nozzle ................65 5.3 Turbulent Cavitating Flow over a ClarkY Hydrofoil.......................................73 6 CONCLUSIONS AND FUTURE RESEACH .....................................................92 6.1 Conclusions of Present Research...................................... ........................ 92 6.2 Future R research D irections........................................................ ............... 94 LIST OF REFEREN CES ............................................................................. 95 BIOGRAPHICAL SKETCH ............................................................. ............... 103 LIST OF TABLES Table page 21 Summary of the hybrid RANS/LES studies ........................................................ 15 22 Parameters used in the computation................................ ................................ 19 23 Comparisons of different turbulence models ................................. ............... 24 31 Overview of the compressible cavitation studies................. ............... ...............37 32 Overview of cavitation on pumps, propellers, inducers and turbine blades ............40 51 Timeaveraged liquid volume fraction v/s pressuredensity correlation at multiple points inside the cavity, original IDM with SoS1, LSM.......................... 63 52 Timeaveraged liquid volume fraction v/s pressuredensity correlation at multiple points inside the cavity, original IDM with SoS2, LSM.......................... 63 53 Comparisons of Strouhal number of original and modified IDM with SoS2A ......72 54 Timeaveraged cavity leading and trailing positions of different turbulence models, modified IDM with SoS2A, AOA=5 degrees ............................ 78 55 Timeaveraged cavity leading and trailing positions of different turbulence models, modified IDM with SoS2A, AoA=8 degrees .............................78 56 Comparison of mean CL and CD, LSM, Cloud cavitation c = 0.80, AoA=8 d e g re e s ...................................... ................................... ................ 9 1 LIST OF FIGURES Figure p 11 Different types of cavitation visualization ...................................... ............... 3 12 Isothermal harmonic speedofsound in the twophase mixture .............................6 21 Geom etry configuration of square cylinder ......................................... ................ 19 22 Streamline snapshot on fine grid with A= 0.15D ........................................20 23 Pressure behaviors at the reference point on the fine grid ....................................21 24 Timeaveraged Uvelocity along the horizontal centerline behind the cylinder......23 25 Snapshots of velocity contour. Color raster plot of axial velocities .....................23 26 Streamline snapshots of two turbulence models on fine grid ..............................24 27 Timeaveraged Uvelocity along y at x/D=0.0................................. ... ................ 26 28 Timeaveraged velocities along y at x/D=1.0 .................................................. 26 29 Tim eaveraged uvelocity along y direction ................................. ................26 210 M ean kinetic energy on different grids ........................................ .....................28 211 Comparisons of different filter sizes on timeaveraged vvelocity along y at x/D = 1.0 .............. .................................... ................ ......... 29 212 Comparisons of different filter sizes on kinetic energy .......................................29 213 Timeaveraged viscosity contours of different turbulence models on fine grid.......30 31 Sketch of a cavity in homogeneous flow ...................................... ............... 42 32 Interface vector sketch in a CV ...................................................... .............. 43 51 Valve geom etry ................................. ... ..... ........... ......... 54 52 Computational domains and boundary conditions................................................55 53 Density contour lines of the steady state solution, original IDM with SoS1, L S M ................................................... ...................... ................ 5 6 54 Middle section density contours at different time instants, original IDM with SoS1, LSM .............................................................................. 57 55 Time evolutions at different locations, original IDM with SoS1, LSM .................58 56 Projected 2D streamlines at middle plane and experimental observation ..............59 57 Middle section density contours of different SoS at different instant time, original ID M L SM ............................................ ......................... 61 58 Two different SoS in hollowjet valve flow, original IDM, LSM ...........................63 59 Pressure time evolutions of different SoS, original IDM, LSM.............................. 64 510 Time evolution and spectrum of pressure and density of two SoS definitions at a point at the cavitation vicinity, original IDM LSM ............................................. 64 511 Timeaveraged eddy viscosity contours of different grids, original IDM with S o S 2 c = 1.9 8 .......................................................................67 512 Timeaveraged vapor volume fraction comparisons of different grids, original ID M w ith S o S 2 c = 1.9 8 ........................................................... .....................6 8 513 Timeaveraged comparisons of different turbulence models on fine grid with A = 0.25Lc,, original IDM with SoS2, c =1.98 ................. ............... 68 514 Instantaneous profiles on fine grid with A = 0.25Lc original IDM with SoS2, U = 1.9 8 ............................................................................. 6 9 515 Pressure contours and streamlines comparison of two turbulence models on fine grid, original IDM with SoS2, c = 1.98 ...................................... ............... 70 516 Pressure evolutions of original and modified IDM with SoS2A at a reference p o in t ..............................................................................7 2 517 Cavity shape and recirculation zone during cycling of original and modified ID M ......... ............................... ... ...................................................7 2 518 Timeaveraged volume fraction and velocity comparisons of original and m o dified ID M .................................................................... .. 7 3 519 ClarkY geometry sketch and Grid blocks.........................................................76 520 Grid sensitivity of the timeaveraged u and vvelocity, LSM, AoA=5 .................77 521 Timeaveraged volume fraction contours and streamlines of different turbulence m odels, A oA = 5 ................................................................ ................. 7 8 522 Timeaveraged volume fraction contours and streamlines of different turbulence models, AoA=8 ........................... .................. .......... .............79 523 Time evolutions of cloud cavitation a = 0.55, AoA=50.......................................80 524 Time evolutions of cloud cavitation a = 0.80, AoA=8 .................................. 81 525 Cavity stage comparisons, cloud cavitation a = 0.80, AoA=8 ..............................83 526 Timeaveraged u and vvelocities of two turbulence models, AoA=5 ..................84 527 Timeaveraged u and vvelocities of two turbulence models, AoA=80 ..................86 528 Timeaveraged lift and drag coefficients comparisons ........................................88 529 Timeaveraged viscosity contours, AoA=50 ................................ .............. 90 530 Timeaveraged uvelocity of different SoS treatments, cloud cavitation c = 0.80, A oA =8 .............. .............................................. ..... .............. 91 c C C C1 C'2 Cp Cdest Cprod k m m P q U u U t T xi y yp LIST OF SYMBOLS speedofsound arbitrary 0(1) constant pressure coefficient turbulence model constants empirical constant in the evaporation term empirical constant in the condensation term turbulent kinetic energy evaporation rate condensation rate pressure source term velocity components in Cartesian coordinates nondimensional velocity magnitude of the horizontal component of velocity time, mean flow time scale temperature Cartesian coordinates nondimensional normal distance from the wall normal distance of the first grid point away from the wall a volume fraction 5ij Kronecker delta function F turbulent dissipation rate t laminar viscosity Pt turbulent viscosity , r, curvilinear coordinates generalized dependent variable r Reynolds stress r, wall shear stress vt kinematic viscosity Pm mixture density a cavitation parameter ok, 7 turbulence model constants Subscripts, Superscripts I interface L, 1 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 oo free stream * predicted value 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 FILTERBASED MODELING OF UNSTEADY TURBULENT CAVITATING FLOW COMPUTATIONS By Jiongyang Wu August 2005 Chair: Wei Shyy Major Department: Mechanical and Aerospace Engineering Cavitation plays an important role in the design and operation of fluid machinery and devices because it causes performance degradation, noise, vibration, and erosion. Cavitation involves complex phasechange dynamics, large density ratio between phases, and multiple time scales. Noticeable achievements have been made in employing homogeneous twophase NavierStokes equations for cavitating computations in computational and modeling strategies. However, these issues pose challenges with respect to accuracy, stability, efficiency and robustness because of the complex unsteady interaction associated with cavitation dynamics and turbulence. The present study focuses on developing and assessing computational modeling techniques to provide better insight into unsteady turbulent cavitation dynamics. The ensembleaveraged NavierStokes equations, along with a volume fraction transport equation for cavitation and turbulence closure, are employed. To ensure stability and convergence with good efficiency and accuracy, the pressurebased Pressure Implicit Splitting of Operators (PISO) algorithm is adopted for timedependent computations. The merits of recent transport equationbased cavitation models are first reexamined. To account for the liquidvapor mixture compressibility, different numerical approximations of speedofsound are further investigated and generalized. To enhance the generality and capability of the recent interfacial dynamicsbased cavitation model (IDM), we present an improved approximation for the interfacial velocity via phase transformation. In turbulence modeling, a filterbased model (FBM) derived from the k e twoequation model, an easily deployable conditional averaging method aimed at improving unsteady simulations, is introduced. The cavitation and turbulence models are assessed by unsteady simulations in various geometries including square cylinder, convergentdivergent nozzle, ClarkY hydrofoil, and hollowjet valve. The FBM reduces eddy viscosity and captures better unsteady features in singlephase flow, and yields stronger timedependency in cavitating flows, than the original k model. Various cavitation models show comparable steady state pressure distributions but exhibit substantial variations in unsteady computations. The influence of speedofsound treatments on the outcome of unsteady cavitating flows is documented. By assessing lift and drag coefficients, pressure and velocity distributions from inception to cloud cavitation regimes, the present approach can predict the major flow features with reasonable agreement to experimental data. CHAPTER 1 INTRODUCTION 1.1 Background of Cavitation In liquid flow, if the pressure drops below the vapor pressure, the liquid is unable to withstand the tensile stress and then adjusts its thermodynamic state by forming vapor filled twophase mixed cavities. This phenomenon is known as cavitation (Batchelor 1967). It can occur in a wide variety of fluid machinery components, such as nozzles, injectors, marine propellers and hydrofoils (Knapp et al. 1970). Cavitation is commonly associated with undesired effects, such as pressure fluctuation, noise, vibration and erosion, which can lead to performance reduction, and can damage the hydrodynamic surface. Apart from these negative effects, one known benefit of cavitation is that it can reduce friction drag (Lecoffre 1999, Wang et al. 2001). For example, in highspeed hydrofoil boats, supercavitating propellers or vehicles the gaseous cavity enveloping the external body surface can provide a shield isolated from the liquid, which helps to cut down the friction. 1.1.1 Cavitation Types in Fluids Different kinds of cavitation can be observed depending on flow conditions, fluid properties and different geometries. Each kind of cavitation has characteristics that distinguish it from the others. Major types of cavitation are briefly described below: Traveling cavitation: Individual transient cavities or bubbles form in the liquid, expand or shrink, and then move downstream (Knapp et al. 1970). Typically, traveling cavitation is observed on hydrofoils at a small angle of attack. The number of nuclei present in the upcoming flow highly affects the geometries of the bubbles (Lecoffre 1999) (Figure 11A). Cloud cavitation: Cloud cavitation is periodically caused by vorticity shed into the flow field. It can associate with strong vibration, noise and erosion (Kawanami et al. 1997) (Figure 11B). Sheet cavitation: Sheet cavitation is a fixed, attached cavity or pocket cavitation in a quasisteady sense (Knapp et al. 1970) (Figure 11C). The interface between the liquid and the vapor varies with flows. Supercavitation: In Supercavitation the cavity appears as the envelope of the whole solid body, and can be observed when underwater vehicles operate at very high speeds or in projectiles with a speed of 500 m/s to 1500 m/s (Kirschner 2001). A typical supercavitating hydrofoil is shown in Figure 11D. Vortex cavitation: Vortex cavities form in the cores of vortices in regions of high shear. They can occur on the tips of rotating blades, and can also appear in the separation zones of bluff bodies (Knapp et al. 1970) (Figure 11E). 1.1.2 Cavitation Inception and Parameter The criterion for cavitation inception based on a static approach can be formulated as Eq. (1.1), with P as local pressure and Pv as vapor pressure P A parameter is needed to describe the flow condition relative to those for cavitating flows and to obtain a unique value for each set of dynamically similar cavitation conditions. In cavitation terminology a denotes the parameter of similitude a= (1.2) p2 /2 In Eq. (1.2), P, V, are the reference pressure and velocity, respectively. Usually, they take the infinity values. Unfortunately, a is not a stringent parameter of similitude for cavitating flows. It is a necessary but not a sufficient condition. One reason is that nuclei in free stream flow can impact the twophase flow structures (Lecoffre 1999). Figure 11. Different types of cavitation visualization (Franc et al. 1995). A) Traveling cavitation. B) Cloud cavitation. C) Sheet Cavitation. D) Supercavitation. E) Vortex cavitation 1.2 Research Motivations and Objectives 1.2.1 Challenges and Motivations Cavitation occurs in various engineering systems, such as pumps, hydrofoils and underwater bodies. It is especially prone to occur when fluid machinery operates at off design conditions. Cavitation may cause not only the degradation of the machine performance but flow instability, noise, vibration, and surface damage. Phenomenologically, cavitation often involves complex interactions of turbulence and phasechange dynamics, large density variation between phases up to a ratio of 1000, multiple time scales, and pressure fluctuations. These physical mechanisms are not well understood because of the complex unsteady flow structures associated with cavitation dynamics and turbulence. There are significant computational issues regarding to accuracy, stability, efficiency and robustness of numerical algorithms for turbulent unsteady cavitating flows. Therefore, turbulent cavitating flow computations need to address both turbulence and cavitation modeling issues. Noticeable efforts have been made in employing Navier Stokes equations for simulations in this field. Among the various modeling approaches, the transport equationbased cavitation models (TEM) have received more interest, and both steady and unsteady flow computations have been reported (Singhal et al. 1997, Merkle et al. 1998, Kunz et al. 2000, Ahuja et al. 2001, Venkateswaran et al. 2002, Senocak and Shyy 2002a,b, 2003, 2004a,b). Senocak and Shyy (2002b, 2004a) developed an interfacial dynamicsbased cavitation model (IDM) accounting for cavitation dynamics. Also, Senocak and Shyy (2002b, 2004a) and Wu et al. (2003b) assessed the merits of alternative TEM models. They showed that for steadystate computations, various cavitation models produce comparable pressure distributions. Despite the good agreement, noticeable differences have been observed in the predicted density field, especially in timedependent computations. In IDM, an empirical factor is used to construct the interfacial velocity in timedependent computations. Then the cavity interfacial velocity is linked to the local fluid velocity. Such an approach lacks generality because the interfacial velocity is supposed to be a function of the phasechange process. A much closer investigation is needed to further understand differences in model performance, especially in unsteady cavitating flows. The differences mentioned above also imply that the compressibility characteristics embodied in each cavitation model are different. This aspect can be significant in the unsteady flow computations because the speedofsound directly and substantially affects the cavity timedependent features (Senocak and Shyy 2003, 2004b). Actually, the local speedofsound of the twophase mixture can reduce by an order of magnitude from either value of pure liquid or vapor. It then leads to large enough Mach number variations around the cavity to make the compressibility effect substantial. The harmonic expression for speedofsound in an isothermal twophase mixture (Venkateswaran et al. 2002, Ahuja et al. 2003) can be presented as 1 a 1a C= P( g + g) (1.3) m PgCg P1CI where p,, Pg, p, are mixture density, vapor density, liquid density, respectively, ag is the vapor fraction, and Cm, Cg, C1 are speedofsound in mixture, in pure vapor, and in pure liquid, separately. The behavior of Cm is plotted in Figure 12 as a function of the vapor volume fraction. The speedofsound drops dramatically over a wide range of the mixture regime. Even though the bulk flow is of a very low Mach number flow, the local Mach number in the interface region may become large. Because of the lack of a dependable equation of state for the mixture sound propagation, Senocak and Shyy (2003) presented two different numerical treatments of speedofsound, namely SoS1 and SoS2, and found that they had an impact on the unsteady cavitating behavior in the convergentdivergent nozzles. They also suggested that SoS2 was more likely to produce the correct unsteady behaviors. This underscores that careful and accurate handling of the speedofsound is important. Currently, a robust compressible model for cavitating flows is still the subject of research. Pure liquid 1400 1200 0 t 1000 CO 800 g 600 400 Pure vapor 200 Flow in this range can be transonic even supersonic 20 Flow, i i I I iI 0 025 05 075 1 Vapor volume fraction Figure 12. Isothermal harmonic speedofsound in the twophase mixture Besides cavitation modeling, the turbulence model can significantly influence the cavitating flow structures. Serious implications of turbulence modeling on cavitating flows were recently revealed by researchers (Senocak and Shyy 2002a, Kunz et al. 2003, CoutierDelgosha et al. 2003, Wu et al. 2003b). They reported that high eddy viscosity of the original LaunderSpalding version of the k Reynoldsaveraged NavierStokes (RANS) model (Launder and Spalding 1974) can dampen the vortex shedding motion and excessively attenuate the cavitation instabilities. Wu et al. (2003a) conducted the nonequilibrium modification, including stationary and nonstationary, in the cavitating flows over a hollowjet valve, and observed no striking impact on the incipient cavitation. Consequently, simulation of phenomena such as periodic cavity inception and detachment requires improved modeling approaches. The Large Eddy Simulation (LES) approach, originally proposed by Smagorinsky (1963) and refined by many researchers (Piomelli 1999, Moin 2002, Sagaut 2003) is an actively pursued route to simulate turbulent flows. However, it is fundamentally difficult to find a gridindependent LES solution unless one explicitly assigns a filter scale (Moin 2002), making the stateofthe art immature for cavitating flow computations. On the other hand, attempts have also been made to employ the information obtained from Direct Numerical Simulations (DNS) to supplement lower order models (e.g. Sandham et al. 2001). To our knowledge, no efforts have been reported to employ LES or DNS for turbulent cavitating flows of practical interest. Recently, efforts have been made to combine the filter concept and the RANS model in singlephase (Mavridis et al. 1998, Batten and Goldberg 2002, Nichols and Nelson 2003, Nakayama 2002, Breuer et al. 2003) and, recently and even more preliminary, cavitating flow computations (Kunz et al. 2003). Hence, it is useful to identify ways to improve the predictive capability of the updated cavitation models and the current RANSbased engineering turbulence closures, which can keep the advantages of RANS approaches and can be easily implemented in practical engineering applications with a clear physical concept. This can help better capture unsteady characters (e.g. the shedding hardly achieved in current turbulent cavitating flow simulations), and can provide better insight into the interactions between cavitation and turbulence models. This can also help establish a dependable, robust and accurate computational CFD tool to analyze and minimize the cavitation effects in the design stage of fluid machinery components or systems. 1.2.2 Research Objectives Our goal is to develop and assess computational modeling techniques to provide better insight into unsteady turbulent cavitation dynamics, aimed at improving the handling of the abovementioned issues in practical engineering problems. Objective 1: Investigate alternative turbulence modeling strategies in the context of the Favreaveraged NavierStokes approach by doing the following: Develop a filterbased turbulence model (FBM) as the alternative to the original LaunderSpalding Model, both employing the k e twoequation closure. Use the FBM in the cavitation simulations to provide insight into the dynamic interactions between turbulence and cavitation. Objective 2: Enhance the cavitation model prediction capability based on the recent transport equationbased model (TEM) by doing the following: Further examine the merits of the recent achievements in the TEM, focusing on the interplay of the twophase compressibility and turbulence by using correlation and spectral analysis with the direct evaluation by experimental measurements, and generalize the speedofsound numerical treatment. Adopt an improved approximation to construct the interfacial velocity by accounting for phase transformation based on the recently developed interfacial dynamicsbased cavitation model (IDM). CHAPTER 2 NAVIERSTOKES EQUATIONS AND TURBULENCE MODELS 2.1 Governing Equations The NavierStokes equations in their conservative form governing a Newtonian fluid without body forces and heat transfer are presented below in the Cartesian coordinates ap 0 puj) + = 0 (2.1) at 9x (pu')+ (puyuj)= OPK+ [ + a 2 0u, (2.2) at ax 9x, 9x 9x ax 3 Ox, The viscous stress tensor is given by S= P + 2 (2.3) ay Ox x, 3 ax, In theory, direct numerical simulation (DNS) can be used to solve Equations (2.1) (2.3) (Rogallo and Moin 1984). However, in practice the limited computing resources prevent us from pursuing such endeavors for highReynoldsnumber flows. Since we are most interested in predicting the mean properties of the turbulent flow, we can first conduct an averaging procedure to simplify the content of the equations. For time dependent flow computations, ensemble averaging is an appropriate conceptual framework. Specifically, to avoid the additional terms involving the products of fluctuations between density and other variables in variable density flows in Reynolds time averaging, Favreaveraging is preferred (Favre 1965) in the following form. Further details can be obtained from many references (e.g. Shyy et al. 1997, Wilcox 1998) = ; p = = +0"; pf"= 0 (2.4) p Then the Favreaveraged NS equations become + V. (p,) = 0 (2.5) at + V. (puu"i)= VP + V. (rz pu") (2.6) 8u, Oui 2 f/ u" B u" 2 c" r, = P( 2 + ' 3,)+Jp( 2&t+~ ) (2.7) adx ax 3 av ax ax 3 a Note that the viscosity fluctuation is neglected. The nonlinear terms (pu"u" ), namely the Reynolds stresses, need additional modeling. The Boussinesq's eddy viscosity hypothesis for turbulence closure leads to the following form j, = =a t( + )39, ( + pk) (2.8) ax 9xI 3 ay Finally, the Favreaveraged NS equations reach the following form S+ V. (p ) = 0 (2.9) at P+V. (p ) = VP+V.(, +r+R) (2.10) R fi7 O 2ai j, + r, = (U +U,)() (2.11) cx cx, 3 c Compared with the original NS equations, Favreaveraged NS equations have the same apparent structure. The above equations can be cast in generalized curvilinear coordinates. The procedure for the transformation of these governing equations is well established (Shyy 1994, Thakur et al. 2002). 2.2 Turbulence Model There are several types of turbulence models. Commonly the classification can be presented as: (a) Algebraic (zeroequation) models (Baldwin and Lomax 1978); (b) One equation models (Bradshaw et al. 1967); (c) Twoequation model, such as k two equation (Launder and Spalding 1974), k o twoequation (Saffman 1970); (d) Second moment closure models (Launder et al. 1975). Among the above models, the kg turbulence model (Launder and Spalding 1974) has been popular because it is computationally tractable with deficiencies reasonably well documented (Shyy et al. 1997, Wilcox 1998). In this model, two partial differential equations accounting for the transport of turbulent kinetic energy k and for dissipation rate s are solved. The following transport equations follow the concept of Launder and Spalding (1974) and are commonly adopted as 1 Au" pk = pu"u; pE = (2.12) a(pk) a(pi k) a a .) t+ = P [( + + ( (2.13) tj xj x ok 8j (pe) d(P ) E E2 A a+ =C P, Clp + [ ) ] (2.14) at Qx k k x Ua Oxj where ok and os are the turbulence model constants, Ce; and Cs2 are the turbulence model parameters that regulate the production and destruction of dissipation of turbulence kinetic energy, respectively. The turbulent production term P, and the turbulent viscosity are defined as p R ; UtP (2p15) Wall functions are used to address the effect of wall boundaries along with the k turbulence model. The empirical coefficients originally proposed by Launder and Spalding (1974), assuming local equilibrium between production and dissipation of turbulent kinetic energy, as given below C,, =1.44, C,, =1.92, = 1.3, k = 1.0 (2.16) In the following, we call it the LaunderSpalding model (LSM). 2.3 FilterBased RANS Model Turbulence plays a very important role in flow phenomena, especially since the Reynolds numbers are high in practical engineering problems. The Reynoldsaveraged NavierStokes (RANS) and, for variable density flows, the corresponding Favreaveraged NavierStokes (FANS) models, such as the ks twoequation closure, have been very popular in providing good prediction for a wide variety of flows with presently available computational resources. However, RANS models describe flows in a statistical sense typically leading to timeaveraged pressure and velocity fields. Generally these approaches are not able to distinguish between quasiperiodic largescale and turbulent chaotic smallscale features of the flow field. The representation may lose the unsteady characteristics when the flow field is governed by both phenomena, even with the help of the nonequilibrium modifications on the set of empirical constants (Wu et al. 2003a). It is clear that the statistical turbulence models have difficulties with the complex phenomena, such as flows past bluff bodies which involve separation and reattachment, unsteady vortex shedding and bimodal behaviors, high turbulence, largescale turbulent structures as well as curved shear layers (Frank and Rodi 1993, Rodi 1997, Shyy et al. 1997). On the other hand, Large Eddy Simulation (LES) operates with unsteady fields of physical values. Spatial filtering is applied instead of averaging in time or ensemble, and turbulent stresses are divided into resolved and modeled parts, such as subgridscale models (SGS) with Smagorinsky's hypothesis (Breuer et al. 2003). Only the large energycontaining eddies are numerically resolved, accomplished by filtering out the high frequency component of the flow fields and using the lowpassfiltered form of the NS equations to solve for the largescale component (Kosovic 1997). In recent years, attempts have been made to remedy the gap between RANS and LES, called the Hybrid RANS/LES methods (Nakayma and Vengadesan 2001, Breuer et al. 2003, Nelson and Nichols 2003). Various strategies have been investigated as summarized below Adhoc models (Kato and Launder 1993, Bosch and Rodi 1996,1998): using a rotation parameter to replace the original production term. Filter RANS/LES model (Koutmos and Mavirdis 1997): combining elements from both LES and standard eddyviscosity approaches, by comparing the characteristic length with the mesh sizespatial filter to reconstruct the viscosity. Multiple Timescale (MTS) Method (Hanjalic et al. 1980, Nichols and Nelson 2003): dividing the turbulent energy spectrum as two parts and breaking the standard ks equations into two sets of equations. Detached Eddy Simulation (DES) (Spalart 1997, Roy et al. 2003): keeping the whole boundary layer (attached eddies) to a RANS model and only the separated regions (detached eddies) to LES. 2.3.1 Literature Review Firstly, we will review the recent studies related to the hybrid RANS/LES studies according to the above categories with a summary in Table 21. Bosch and Rodi (1996) applied the adhoc model, which used a rotation term to reduce the turbulent production, to simulate the vortex shedding past a square cylinder near a wall. Compared with the standard turbulence model, this modification produced better unsteady behavior and obtained good agreement with the experimental data. Following their previous study, Bosch and Rodi (1998) adopted a 2D ensembleaveraged unsteady NavierStokes equation, with the adhoc model to compute the vortex shedding past a square cylinder. The numerical results agreed well with experimental measurement and other similar numerical results. They also carried out other versions of turbulence models in the same configurations for comparisons. They concluded the present modification was better. Rodi (1997) simulated turbulent flows over two basic bluff bodies, 2D square cylinder and 3D surfacemountedcube, using different Reynolds numbers and a variety of LES and RANS methods. The various calculation results generally agreed with detailed experimental data. Assessment was given based on performance, cost and the potential of the various methods based on the comparison with the measurement. Koutmos and Mavridis (1997) combined LES and the standard kg models to formulate the eddyviscosity by comparing a mesh size A with the characteristic length in the flow field. Then the turbulent viscosity was constructed in two different ways. The calculation of unsteady separated flows of square cylinder wake and backwardfacing step recirculating flows under low (Re=5000) and highReynolds number (Re=37000) conditions agreed well with the experimental measurements. Nagano et al. (1997) developed a low Reynolds number multiple timescale turbulence model (MTS) which separated the turbulent energy spectrum into two parts, namely, production and transfer according to the wavenumber. Then the eddyviscosity was approximated by solving both wall and homogeneous shear flows. The test results compared well with the DNS and experimental data. They concluded for homogeneous shear flow the difference between the new model and the standard ks model from the estimation of the characteristic timescales, not from the discrepancy in the eddy approximation. Nichols and Nelson (2003) employed different turbulence models, including RANS, MTS and DES to simulate several unsteady flows. Based on the comparison with experimental data, they made the assessment of different models, and suggested that the MTS hybrid RANS/LES model needed more investigations in term of grid and timestep sensitivities. Roy et al. (2003) examined Detached Eddy Simulation (DES) and RANS turbulence modeling in incompressible flow over a square cylinder. They found that the 2D and 3D simulations using DES are almost identical, and the results also compared well with experimental data, while the steadystate RANS significantly overpredicted the recirculating vortex behind the cylinder. Table 21. Summary of the hybrid RANS/LES studies Author and year Methods Conclusions Bosch and Rodi Ensembleaveraged NS with adhoc Capture better unsteady 1996, 1998 modification in turbulence closure shedding than RANS and agreed well with experimental data Table 21. Continued Rodi 1997 RANS with LES as turbulence closure Generally agreed with experiment measurements and found 2D could not resolve the 3D effect Koutmos and Ensembleaveraged NS with a spatial Resolved good unsteady Mavridis 1997 filtering form LES features in square cylinder and backward facing step under low and high Reynolds numbers Nagano et al. A lowReynolds number MTS Agreed with the DNS and 1997 turbulence model in homogeneous shear experiment data flow Nichols and RANS, Hybrid RANS/LES, DES Hybrid models is better than Nelson 2003 RANS with comparable accuracy with LES Roy et al. 2003 NS equation with a DES turbulence Compared well with closure experimental data 2.3.2 FilterBased Model Concept (FBM) In the RANS models developed for steady state flows, the turbulent length scales predicted by the model extend over a large part of the flow domain. By imposing a filter on the flow, the turbulent scales smaller than the filter will not be resolved. When the filter size is smaller than the length scales returned by the RANS models, this will allow the development of flow structures that are not dissipated by the modeled effective viscosity. The subfilter flow may be characterized by transport equations for turbulent energy, dissipation, and Reynolds stresses. In the present example, we choose to apply the LSM (Launder and Spalding 1974) as the corresponding RANS model. The filtering operation will be controlled by the size of the imposed filter A and the size of the RANS length scale lINs. More detailed can be obtained from Johansen et al. (2004). The brief concept is given below. We start from the RANS length scale 1/los =c3 ki (2.17) We may now construct filtered eddy viscosity in the following general form V C412 C 2 Akc R4NS A. vt =C =C NSf(C3 3) =3 V (C ) (2.18) where lf is the turbulent length scale survived during the filtering operation. The scaling function f must impose the filter, and have limiting properties such as f(C3 ) = 1 exp(C3 (2.19) k3/2 k3/2 Here k and E are the nonresolved turbulent energy and dissipation rate separately, corresponding to unresolved turbulent length scales ,ff 1, = 1 exp(C3 k ) (2.20) The proposed model becomes identical to the RANS model in the extremely coarse filter. In the case of a fine filter, the turbulent length scale is controlled only by the imposed filter, and the LES type of model is obtained under this condition. The model is expected to have the following properties If filter size A is identical to cell size and grid Reynolds numbers are sufficiently small (Kolmogorov scales) the DNS model is recovered. If filter size becomes large, the RANS model is recovered. The statistical understanding of the advancement of the averaged flow during one timestep implies that the time step itself should be a part of the model. The filter should be almost independent to the grid in a specific computation, though it may vary with different geometries (different characteristic lengths). 2.3.3 Modeling Implementation The k E twoequation is adopted in similar formula as Eq.(2.13)(2.14). In Eq.(2.15), the turbulent viscosity modeled with a filter by the scaling function of Eq.(2.19) leads to the following Ac k2 = 0.09.Min[1,c3 cA ] (2.21) k3/2 E with C3 1, here we choose C3 = 1.0. This choice helped to assure that in near wall A8g nodes the scaling function f = Min[l,c3 k ] will always return f = 1.0. Then the use of standard wall functions is fully justified as the standard k e model. In the limit of a very large filter size A, the viscosity becomes v, = 0.09 k2 / e and the standard k e model is recovered. In the limit of a filter much smaller than the turbulent scale k3/2 I/, the viscosity model becomes v, = 0.09 C3Ak1/2 (2.22) Then the FBM is identical to the oneequation LES models of Schumann (1975) and Yoshizawa (1993). 2.4 Assessing FBM in TimeDependent SinglePhase Flow For the filterbased model (FBM), we first carry out the studies on the singlephase simulations of vortex shedding past a square cylinder, and compare with experimental measurements (Lyn and Rodi 1994, Lyn et al. 1995). The parameters for the computations are given in Table 22. The square cylinder, with a height D, is located at the center of a water channel. Figure 21 shows the geometry of the experimental setup of Lyn and Rodi (1994) and Lyn et al. (1995), which is selected to guide the investigation, and the nondimensionalized coordinates. It consists of a 2D square cylinder inside a channel. To reduce numerical errors in the vicinity of the cylinder, we adopt uniform grid spacing within the 4Dx3D domain surrounding the cylinder. Outside this block, the grid is slightly expanded toward the edges of the 19 computational domain. All variables are nondimensionalized by the free stream velocity and the cylinder height. The fluid properties are held unchanged and the Mach number is zero. At the inlet, the mean velocity U,, is uniform and follows the horizontal direction. The inlet turbulent intensity is 2%. Based on the definition of the turbulent viscosity formula adopted in the LSM and by assigning the turbulent Reynolds number to be 200, we determine the inlet dissipation rate. The flow variables are extrapolated at the outlet. The wall function (Shyy et al. 1997) is employed for the solid boundary treatment. Table 22. Parameters used in the computation U.n (m/s) p (kg/m3) p (kg/m.s) Re D (m) k (m2/s2) 2 (mZ/s3) 0.535 1000 1.002x103 21357 0.04 1.7174x10 4 2.5x10 6 Solid boundary D=0.04 020 0.56 6.5D 1.0 4  I 3c Outlet Outlet 1 I Ceterline 0 L Soihudr 1 2 31 4 1 5 1 6 1 Solid boundary 0 1 2 3 4 5 6 x/D A B Figure 21. Geometry configuration of square cylinder. A) Computational domain. B) Nondimensional coordinates To investigate the effect of grid resolution on numerical accuracy, we used three levels of grids including fine, intermediate and coarse grids, which have 25, 20 and 10 nodes on each side of the cylinder respectively. For comparison, the LaunderSpalding model (LSM) is also carried out on coarse and fine grids. To investigate the sensitivity, we use four different filter sizes: 0.15D, 0.3D, 0.6D and 0.9D on coarse grid (10 intervals), and 0.15D, 0.25D, 0.6D and 0.9D on fine grid (25 intervals). Unless explicitly mentioned, A= 0.15D is used as the reference filter size. Inlet If (A E) / k3/2 exceeds 1.0, the filter scaling function f Min[1.0, (A )/ k3/2] will take a value of 1.0. This treatment enables one to apply the wall function of LSM for the solid wall treatment, which is confirmed by the outcome presented in Figure 22. A pressure reference point was located at position x/D=0.0 and y/D=0.50, which is the midpoint of the cylinder upper wall. For the fine grid we show predicted pressure development at the wall reference point in Figure 23A. We see a modulated and inexact periodic signal that qualitatively agrees with the experiments (Lyn and Rodi, 1994). The FBM produces pressure oscillations corresponding to amplitudes of pressure coefficient of approximately 2.5, while pressure coefficient amplitude predicted by the LSM varies from approximately 0.1 to less than 0.05. We also found that the LSM tends to be more time independent on the finer grid and the oscillations die out eventually. The FFT of the Cp by filterbased model (Figure 23A) shows multiple frequencies with a dominant one at around 2.15 (Figure 23B). A B Figure 22. Streamline snapshot on fine grid with A = 0.15D Red shadow indicates f = Min(.0, 3 ) = 1.0 for recovering the LSM. At other regions the filter function is employed. A) Zoom in. B) Full view SLaunderSpalding 1800 05 Filterbased 1600 1 1400 5 11200 C 1000 2 C 800 2 5 600 400 200 5 10 15 20 1 2 3 4 5 t s AFrequency A B Figure 23. Pressure behaviors at the reference point on the fine grid. A) Pressure coefficients Cp time evolutions, LSM and FBM. B) Cp FFT, FBM Figure 24 shows the timeaveraged horizontal velocities along the centerline behind the cylinder. The LSM predicts too long wake lengths 4 behind the cylinder at about 3.0D for both grids (coarse and fine grids), which are almost identical to the numerical results by Rodi (1997). Further results of the FBM, using a constant filter size of A = 0.15D are shown in the same figure. For the intermediate and fine grids, the FBM results quantitatively agree with the experimental data of Lyn et al. (1995) including the asymmetric behavior. Even for the coarse grid, the size of the separation zone is well reproduced. However, the reverse velocity in the wake as well as the velocity defect in the remaining part of the wake is slightly underpredicted. In the case of the coarse grid the resolution of the shear layer at the cylinder wall is subcritical and this affects both the onset of each shedding cycle as well as the magnitude of the vorticity in the wake. Flow structures of the solutions at a given time instant on coarse and fine grids, and with LSM and FBM are highlighted in Figure 25. For the FBM, as the grid resolution is refined, the vortex structure becomes more dispersed and less confined in the wake region. In contrast, for the LSM, the impact of the turbulent viscosity is dominant, and the fine grid solution exhibits less fluctuation in time. The streamlines at one time instant for both models are presented in Figure 26. The LSM only catches two separated pockets at the cylinder upper and lower shoulders, and two wakes behind the cylinder. The shedding motion is almost gone. The FBM is able to capture the sharp separation in the shear layer, which agrees well with the experimental observation from Lyn and Rodi (1994) and Lyn et al. (1995). From Table 23 we see that the predicted Strouhal numbers are about 20% higher than the St=0.135 from experiment by Lyn and Rodi (1994). By comparing the present solutions on different grids using the FBM with the same filter size, the variations of the Strouhal number are less than 4%. Hence, the only significant result is the Strouhal number being approximately 20% larger. This overprediction may be due to the 2D effect, compared with the 3D measurement in the experiment. The largescale structures of the flow have a threedimensional nature and we do not expect to reproduce all features of this flow correctly in two dimensions. This will be more pressing as the filter size and grid is reduced and we depart more and more from the LSM. This seems to be consistent with the results from Yu and Kareem (1997) who needed to use a larger Smagorinsky coefficient to reproduce the Strouhal number in 2D compared to their full 3D simulations. Another interesting finding is that the LSM required much smaller time steps, compared with the FBM in order to get stable and convergent solutions using the PISO algorithm. From Table 23 it is seen that the time steps have to be reduced substantially in order to employ the LSM. Hence, the time consumption with the FBM is around half smaller in the present study. 08 / / 06 / /1 '04 02 LSM, coarse LSM, fine S A=0.15D, coarse 0 A  A =0.15D, intermediate 0 2 5 10 15 x/D Figure 24. Timeaveraged Uvelocity along the horizontal centerline behind the cylinder. Experimental data are from Lyn and Rodi (1994) and Lyn et al. (1995) Coarse grid Fine grid a Coar W7 Figure 25. Snapshots of velocity contour. Color raster plot of axial velocities (red is largest, blue lowest). A) FBM: coarse grid, velocity range is 0.311.53. B) FBM: fine grid, velocity range is 0.661.88. C) LSM: coarse grid, velocity range is 0.09 to 1.25. D) LSM: fine grid, velocity range is 0.31 to 1.20 r [00 __C~ ~"x~,~  __I~i: Figure 26. Streamline snapshots of two turbulence models on fine grid. A) LSM. B) FBM Table 23. Comparisons of different turbulence models. /ID is the relative position of the reattachment, measured from the cylinder center (coarse: 10 intervals, intermediate: 20 intervals, fine: 25 intervals on each cylinder face) Model Grid (n, x ny) Filter Size Time Step I/D St dt/(D/Uin) 0.15 D 0.0134 1.22 0.155 0.30 D 0.0134 1.40 0.151 Coarse (162x92).6D 0.0134 2.12 0.143 0.90 D 0.0134 2.73 0.137 FBM Intermediate (290 x190) 0.15 D 0.0134 1.25 0.163 0.15 D 0.0669 1.23 0.161 0.25D 0.0669 1.25 0.148 Fne(30095) 0.60D 0.0669 1.98 0.140 0.90D 0.0669 2.64 0.134 Coarse (162x92) 0.00268 3.03 0.124 LSM Fine (300x195) 0.000803 2.80 0.125 Exp. Lyn et al. (1995)38 01 The predicted velocity profiles are given in Figure 27 to Figure 210. One thing to be noted is that the experiments of Lyn and Rodi (1994) and Lyn et al. (1995) were recorded on a single side, assuming that the ensemble and time averaged flow was symmetrical across the axial symmetry line. As a result, data were recorded only at one side of the symmetry line. Hence, we have mirrored the data to enable the visual I 'i'i \'i i A1 :~fl comparison with the current numerical results. By close inspection of the experimental data, we find that the transversal time averaged velocity on the symmetry line is non vanishing. This indicates that the data are not completely symmetrical. Along the vertical centerline (Figure 27), we see that the LSM is unable to correctly reproduce the separation in the shear layer. The shoulders on the velocity profiles are not captured, presumably due to the high effective viscosity in the incoming flow. The FBM solution on the coarse grid is similar to the fine grid solution, but only as a result of poor resolution of the shear layer. The profiles for the FBM match the experimental data of Lyn and Rodi (1994). The timeaveraged transversal velocity in the wake, one cylinder diameter behind the center of the cylinder, is shown in Figure 28B, where the LSM fails to capture the correct transversal velocity. For both intermediate and fine grids, the FBM again gives very good results. The timeaveraged axial velocity is shown in Figure 28A. Here we see that the FBM solutions become more asymmetrical by grid refinement. The asymmetry in the timeaveraged solutions is clearly seen from Figure 28 to Figure 29, and has been confirmed independently by calculations using a code that employs a staggered grid arrangement. The development of asymmetrical time averaged solutions seemed to be caused by the initial bias of the solution, which is path dependent. Asymmetries in the timeaveraged fields were noted by Sohankar et al. (1999). They reported that the asymmetries became more distinct as the Reynolds number increased. This is consistent with our findings. As suggested by Sohankar et al. (1999), the reason could be the twodimensional geometry that is forced on the flow. In two dimensions wall attachment may be reinforced due to the imposed two dimensionality. This can only be resolved by comparing full 3D computations. 08 _ A=0.15D, coarse  A=0.15D, intermediate 3 06 ' A=0.15D, fine SA Experiment 04 02 I 0 02 04 I 1 0 1 y/D Figure 27. Timeaveraged Uvelocity along y at x/D=0.0. Experimental data are from Lyn and Rodi (1994) and Lyn et al. (1995) I ; A 8 / 'f i" 1 I SI I LSM, coarse / / LSM, fine / A=0.15D, coarse ) / .. A=0.15D, intermediate S A=0.15D,fine \ AA A Experiment 2 1 0 1 2 y/D LSM, coarse A LSM, fine A/ \ A=0.15D, coarse A \ _._._. A=0.15D, intermediate S  A=0.15D, fine S\\ A Experiment A '1^ \ \A \ A 2 1 0 y/D 1 2 Figure 28. Timeaveraged velocities along y at x/D=1.0. Experimental data are from Lyn and Rodi (1994) and Lyn et al. (1995). A) Uvelocity. B) Vvelocity  / A\ \ I' 'I \ A AA I A SLSM, coarse LSM, fine SA=.15D, coarse \ A=0.15D, interme \ A=0.15D,fine Experiment , I I I I I 4 3 2 1 0 1 2 3 4 y/D Figure 29. Timeaveraged uvelocity along y direction. Experimental data are from Lyn and Rodi (1994) and Lyn et al. (1995). A) At x/D=3.0. B) At x/D=8.0 The predicted kinetic energies along the horizontal centerline are compared with experiments in Figure 210. We find that the LSM significantly underpredicts the total kinetic energy, including the periodic and turbulent parts, mainly due to poor resolution of the large scale structures in the wake. The FBM results agree well with the experiment, even in the coarse grid, in terms of the magnitude and the peak position (Figure 210A). The turbulent kinetic energy (the stochastic part only) is underpredicted for both LSM and FBM, as shown in Figure 210B. Compared with the total and turbulent kinetic energy, the turbulent part dominates for the LSM, while the periodic part is dominant for the FBM. Since the FBM overpredicts the turbulent kinetic energy, it slightly over predicts the periodic kinetic energy at the peak. Next, we investigate the effect of the filter size. As mentioned above, we present the results on the coarse grid with four filter sizes 0.15D, 0.3D, 0.6D and 0.9D, and fine grid with four filter size 0.15D, 0.25D, 0.6D and 0.9D. Figure 211 shows the transversal vvelocity of both coarse and fine grids with different filter sizes. The profiles exhibit a clear trend toward the LSM when the filter size increases. Table 23 further demonstrates that the computed Strouhal number and the reattachment length 4/ also move toward the LSM as the filter size increase on the coarse grid. Also, the kinetic energy shows a similar trend with the various filter sizes on both coarse and fine grid (Figure 212). Overall, the solution with A= 0.15D agrees better with the experiment. By inspection of the results it appears that the present filterbased calculations give quite regular solutions. In many calculations of the shedding from a square cylinder, perturbations of the flow are induced by "numerical noise" that may be caused by a large number of different phenomena. Examples of such phenomena are unbounded convective fluxes, reduced numerical order caused by expanded or unstructured grids, and too large time step sizes. In the filterbased model randomness can be added to the solution by applying a random force field, similar to what is used in RenormalizationGroup analyses (Smith and Woodruff, 1998). Without inducing randomness to the flow by inlet conditions or by random forcing we expect that the present model will produce regular oscillating solutions similar to URANS (unsteady RANS) calculations by laccarino et al. (2003). The viscosity contours of LSM and FBM are shown in Figure 213. The LSM predicts a very high viscosity in the upstream of the cylinder, and the distribution is almost symmetric. The FBM significantly reduces the viscosity around 2orders of magnitude. As discussed above, the high effective viscosity in the coming flow can damp out the unsteadiness behind the cylinder. LSM, coarse A LSM,coarse / LSM, fine A LSM, fine 6  A=0.15D, coarse A=0.15D, coarse /  A=0.15D, intermediate A A=0.15D, intermediate 1// . A=0.15D, fine 025 A=0.15D, fine i A Experiment A A Experiment 4 I 02 A 04 F A / \ An M r ,0' A A A A A A 02 / A A 01 A A 02 \ A A "0 05 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 x/D x/D A B Figure 210. Mean kinetic energy on different grids. Experimental data are from Lyn and Rodi (1994) and Lyn et al. (1995). A) Total kinetic energy (mean + turbulent). B) Turbulent kinetic energy (stochastic) LSM, fine A A=0.90D, fine V\  A=0.60D, fine IA \ A=0.25D, fine \/ A=0.15D, fine / \ Experiment V N \ M \\< ,! I\ l 2 1 0 y/D 1 2 A B Figure 211. Comparisons of different filter sizes on timeaveraged vvelocity along y at x/D=1.0. Experimental data are from Lyn and Rodi (1994) and Lyn et al. (1995). A) Coarse grid. B) Fine grid LSM, coarse A=0.90SOD, coarse A=0.60D, coarse .A=0.30D, coarse A:0?0 .....coa I . ....... .. A=0.15D, coar A /" A Experiment \.. A/ ...  ...  / A / ......... V /  V/ 1 2 3 4 5 6 7 x/D  A \ A I '.  / // I I I I 1 2 3 4 x/D LSM, fine A=0.90D, fine A=0.60D, fine  A=0.25D, fine ... A=0.15D, fine A Experiment AA ~'Lr~ 03 025 01 LSM, coarse A=0.SO9D, coarse ..... A=0.60D, coarse  A=0.30D, coarse  A=0.15D, coarse A Experiment A A A   1 2 3 4 x/D 03 025 02 .015 01 5 6 7 5 6 7  LSM, fine  A=0.90D, fine . A=0.60D, fine  A=0.25D, fine ..... A=0.15D, fine A Experiment A A A A A A A A 1 2 3 4 x/D 5 6 7 Figure 212. Comparisons of different filter sizes on kinetic energy. Experimental data are from Lyn and Rodi (1994) and Lyn et al. (1995). A) Coarse grid: total kinetic energy (periodic + turbulent). B) Coarse grid: turbulent kinetic energy (stochastic). C) Fine grid: total kinetic energy (periodic + turbulent). D) Fine grid: turbulent kinetic energy (stochastic) 7" 06 04 021 02 nI U )5  . ie I I        S1.81 .g09 0.182 mo.oos A B Figure 213. Timeaveraged viscosity contours of different turbulence models on fine grid. A) LSM. B) FBM In summary, a filteredNavierStokes model, originated from the LSM, is applied to vortex shedding from a square cylinder. The introduction of the filter led to an effective viscosity that depends on both turbulent quantities and the filter size itself. The method is capable of working with standard wallfunctions, allowing much coarser grids in the boundary layer compared to common LES methods. Presently the use of wall functions is justified as y' values for near wall nodes are greater than 20 using the fine grid (25node). Based on the discussions above, we have the following conclusions Both coarse and fine grids reproduce the time averaged experimental results quantitatively. However, by refining the grid we see improved results for the velocity profiles. For the investigated filter size of A = 0.15D, the solutions on both intermediate and fine grids are in agreement with experimental data, demonstrating that the model can produce better resolutions based on the LSM by allowing the numerical scheme to simulate the fluid physics at the scales where numerical resolutions are satisfactory. The increase of filter size shows that the filterbased model smoothly approaches the LSM. The filterbased model is shown to produce improvement over the LSM for all grids investigated. The Strouhal number of FBM is generally overpredicted. Further investigation is needed to investigate whether it can be resolved by full 3D simulations. 31 Generally, the model is expected to yield better results if full 3D solutions are performed, since the large scale 3D coherent flow structures can be resolved. CHAPTER 3 CAVITATION MODELS In this chapter, the single set of the governing equations for the flow field, including the continuity and momentum equations which were introduced in Chapter2 along with a transport equation for the cavitation, is presented first. Then, the cavitation models utilized for the present study are provided, along with a review of selected recent studies. The above equations coupled with turbulence models, including Launder Spalding model (LSM) and Filterbased model (FBM), presented in the previous chapter will complete the whole system of equations for the turbulent cavitating flow computations. 3.1 Governing Equations The Favreaveraged NavierStokes equations, in their conservative form, are employed for incompressible flows. The cavitation is governed by a volume fraction transport equation. The equations are presented below in the Cartesian coordinates OP. a(Pu, ) + = 0 (3.1) at 9x a(Pmu,) C(Pmu,Uj) op a au, ?uJ 2 al/ + + [(P+,u)( + 5,)] (3.2) at dCx 9x 9x O)x cx 3 (x aa+ (auj) = (ih +i' ) (3.3) at ox where p, is the mixture density, u, is the velocity component in Cartesian coordinates, t is the time, p is the pressure, a, is volume fraction of liquid, ih+ is the condensation rate and li is the evaporation rate, /t is the turbulent viscosity, and / is the laminar viscosity of the mixture defined as /A=/ ,a, +/,(1a,) (3.4) with /u is the laminar viscosity of the liquid and /u is the viscosity of the vapor. The mixture density and turbulent viscosity are respectively defined as below p, = plo,,+ (l a,) (3.5) PmCjik2 =t = (3.6) where p, is the density of the liquid and p, is the density of the vapor, C, is the turbulence model constant and its value is 0.09 generally, k is the turbulence kinetic energy and E is the turbulence dissipation rate. The relevant details of the different turbulence closures have been reviewed in detail in Section 2.22.4. 3.2 Literature Review of Cavitation Studies Cavitation can produce negative effects in fluid machinery components and systems. Details of the existence, extent and effects of cavitation can help to minimize cavitation effects and optimize the designs. Experiments have been conducted in the past few decades for different types of fluid machinery devices and components. Ruggeri and Moor (1969) investigated methods that predicted the performance of pumps under cavitating condition. Typically, the strategy to predict the Net Positive Suction Head (NPSH) was developed. Stutz and Reboud (1997, 2000) studied the twophase flow structure of unsteady sheet cavitation in a convergentdivergent nozzle. Wang (1999) used highspeed camera and Laser Light Sheet (LLS) to observe the cavitation in a hollowjet valve under different openings. Wang et al (2001) studied broad cavitation regimes of turbulent cavitating flows, covering from inception to supercavitation, over a ClarkY hydrofoil under two different angleofattacks. Besides experimentation, using CFD simulation to analyze cavitation phenomena has become convenient and popular with the development of computer hardware and software. A number of cavitation models have been developed. These studies can be put into two categories, namely interface tracking methods and homogeneous equilibrium flow models. A review of the representative studies is presented by Wang et al. (2001). Taking the advantage of homogeneous equilibrium flow theory, the mixture concept is introduced. And a unique set of mass and momentum equations along with turbulence and cavitation models is solved in the whole flow field. Within the homogeneous equilibrium flow theory, two approaches can model the cavitation dynamics. The first one is the arbitrary barotropic equation model, which suggests that the relationship between density and pressure is p = f(p), and the second one is the transport equation based model (TEM). Barotopic equations were proposed by Delannoy and Kueny (1990). They assumed that density was a continuous function of pressure where both pure phases were incompressible, and the phase change could be fitted by a sine curve. They could not correctly produce the unsteady behaviors in the venturi simulation. Arbitrary barotropic equation models (density is only a function of pressure) do not have the potential to capture baroclinic vorticity production because the baroclinic term of the vorticity transport equation yields zero by definition (Senocak 2002, Senocak and Shyy 2004a). In addition to agreement with the experimental study of Gopalan and Katz (2000), the above two references have demonstrated computationally that the baroclinic vorticity generation is important in the closure region. In TEM, a transport equation for either mass or volume fraction, with appropriate source terms to regulate the mass transfer between vapor and liquid phases, is adopted. The 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 cavity bubbles, especially in complex 3D interface situations (Wang et al. 2001). Different modeling concepts embodying qualitatively similar source terms with alternate numerical techniques have been proposed by various researchers (Singhal et al. 1997, Merkle et al. 1998, Kunz et al. 2000, Ahuja et al. 2001, Senocak and Shyy 2003, 2004a,b). Numerically, Singhal et al. (1997, 2002) and Senocak and Shyy (2002a,b, 2003, 2004a,b) utilized pressurebased algorithms, while Merkle et al. (1998) and Kunz et al. (2000) employed the artificial compressibility method. In addition, Vaidyanathan et al. (2003) performed a sensitivity analysis on a transport equationbased cavitation model to optimize the coefficients of its source terms. More recently, attempts have been carried out on numerical simulation of cavitating flows in turbomachines, such as pumps and inducers. Medvitz et al. (2002) utilized the preconditioned twophase NS equations to analyze the performance of cavitating flow in a centrifugal pump. Hosangadi et al. (2004) simulated the cavitating flow in the full geometry of a 3blade Simplex inducer. Couties Delgosha et al. (2005) presented a 3D model for cavitating flow through a 4blade inducer adopting only one bladetoblade passage. However, a complete robust and accurate CFD framework in this field is still a longstanding work. 3.2.1 Cavitation Compressibility Studies Since the flow fields are rich in complexity at the cavity interface region, as pointed out in chapter 1, compressibility effect is one of the major issues in cavitation studies. Nishiyama (1977) developed a linearized subsonic theory for supercavitating hydrofoils to clarify the compressibility correction of Mach number effect in pure water. The 2D & 3D steady characteristics of supercavitating hydrofoils in subsonic flow were compared to those in incompressible flow. The essential differences, including the compressibility effect between vapor and liquid and the corelation between cavitation and Mach number, were shown in detail. Saurel and Cocchi (1999) focused on cavitation in the wake of a highvelocity underwater projectile. They presented a physical model based on the Euler equations in terms of twophase mixture properties. The mathematical closure was achieved by providing state equations for the possible thermodynamic states: compressible liquid, twophase mixture and pure vapor. The model was then solved using a hybrid computational scheme to accurately maintain the property profiles across the discontinuities. The results demonstrated a reasonable agreement compared with the known analytical solutions. Saurel and Lemetayer (2001) proposed a compressible multiphase unconditionally hyperbolic model to deal with a wide range of application: interfaces between compressible materials, shock wave in condensed multiphase mixtures, homogeneous twophase flows and cavitation in liquids. The model did not require a mixture equation of state and was able to provide themodynamic variables for each phase. The results and validations with analytical solutions were provided. Kunz et al. (2000), followed by Venkateswaran et al. (2002), developed a preconditioned timemarching algorithm for the computation of multiphase mixture flows based on carrying out perturbation expansions of the underlying timedependent system. The method was efficient and accurate in both incompressible and compressible flow regimes. However, it was not a fully compressible formulation for the flow fields. It could not be accurate if the bulk flows were supersonic, which indicated that compressibility should be considered even in pure liquid phase. Senocak and Shyy (2004a,b) used the pressurebased method to simulate the cavitating flows through convergentdivergent nozzles. They presented two different treatments of speedofsound (SoS), including SoS1 and SoS2, to build up the relationship between the pressure and density for steady and unsteady cavitating flows. The comparison with the experimental data and other simulation results displayed that SoS2 was more likely better in predicting the unsteady characteristics. Table 31. Overview of the compressible cavitation studies Author and year Methods (Analytic/Numerical) Conclusions(Analytic/Numerical) Nishiyama Proposed linearized theory linearized relationships between 1977 velocity, pressure and sound speed Saurel and A physical model based on the Reasonable agreed with analytic Cocchi 1999 Euler equations, with state solutions and could reliably deal equations with strong shock wave and complex EOS Saurel and An unconditional hyperbolic model, Compared with analytic solutions Lemetayer 2001 with accurate treatment for non And clouded it could deal with a conservative form to ensure mass wide range of application, such as conservation shock wave, homogeneous two phase flow and cavitation. Kunz et al. A preconditioned timemarching Agreed with experiment data 2000, CFD method for isothermal under good precondition, and Venkateswaran multiphase mixture flows, suggested the compressible et al. 2002 associated artificial compressibility treatment could improve the dynamics description than incompressible computations Senocak and Two numerical models for of Agreed well with experimental Shyy 2004a,b speedofsound in steady and data and other numerical unsteady cavitating flows simulations, and demonstrated the capability in unsteady computation. 3.2.2 Cavitation Studies on Fluid Machinery Components and Systems Ruggeri and Moor (1969) investigated similarity methods to predict the performance of pumps under cavitating condition for different temperatures, fluids and operating conditions. Typically, the strategy to predict the NPSH was developed based on two sets of available test data for each pump at the concerned range of operating conditions. Pumps performance under various flow conditions such as discharge coefficient, impellor frequency as assessed for different fluids, such as water, LH2 and butane. Wang (1999) studied the cavitating phenomena in a hollowjet valve under different openings and cavitation conditions, using high speed camera and LLS. Further more, he also studied the induced vibration mechanism, cavitation damage characteristics and ventilation effect in the valve. Athavale and Singhal (2001) presented a homogeneous twophase approach with a transport equation for vapor, and the reduced RayleighPlesset equations for bubble dynamics based on local pressure and velocity. Compared with the experimental measurements, they obtained reasonable predictions of cavitating flows in two typical rocket turbopump elements: inducer and centrifugal impeller. Lee et al. (2001a,b) analyzed the cavitation of the pump inducer sequentially. They found that the forwardswept blade demonstrated more resistance to vortex cavitation than the conventional one. They also observed that the cavity length of surface cavitation at various conditions was closely related with the cavitating number and event duration, but inconsistent with the mode predicted by linear cascade analysis in the crossflow plane at far offdesign points. Medvitz et al. (2002) used the homogeneous twophase RANS equations to analyze the performance of 7blade centrifugal pumps under cavitating conditions. By using the quasi3D analysis, the numerical results were found to be qualitatively comparable with the experimental measurements across a wide range of flow coefficients and cavitation numbers, including offdesign flow, blade cavitation and breakdown. Duttweiler and Brennen (2002) experimentally investigated a previously unrecognized instability on a cavitating propeller in a water tunnel. The cavitation on blades and in the tip vortices was explored through visual observation. The cyclic behavior of the attached blade cavities had strong similarities to that of partial cavity oscillation on single hydrofoils. Furthermore, the reduced frequency of the instability was consistent with the partial cavity instability on a single hydrofoil. Friedrichs and Kosyna (2002) described an experimental investigation of two similar centrifugal pump impellers at low specific speed. The highspeedfilm displayed rotating cavitation over a wide range of part loaded operating points, and illustrated the development of this instability mechanism, which was mainly driven by an iteration of the cavity closure region and the subsequent blades. Hosangadi et al. (2004) numerically simulated the performance characteristics of an inducer using an NS methodology coupled with a twoequation turbulence model. The simulations were performed at a fixed flow rate with different Net suction Specific Speeds (NSS). The numerical results showed the head loss was related to the extent of cavitation blockage. The breakdown NSS and the head loss prediction agreed well with experimental data. The insights provided a sequence of traveling and alternate cavitation phenomena in blade passages. CoutierDelgosha et al. (2005) presented a 3D model for a cavitating flow in 2D venture and in a 4blade inducer with the comparison to experimental data and visuals. They assumed symmetry in the inducer and only considered one bladetoblade channel in the inducer. The quasistatic results showed a consistent agreement with the experiment, but did not catch the performance breakdown. Table 32. Overview of cavitation on pumps, propellers, inducers and turbine blades Author and year Problem study (Experimental) Findings(Experimental) Methods (Numerical) Conclusions(Numerical) Ruggeri and Similarity to predict the pumps Provide strategies to predict the Moore 1969 performance under cavitation pump NPSH conditions Wang 1999 Highspeed camera and LLS in a Cavitation, induced vibration, hollowjet valve damage characteristics and the ventilated cavitation effect Athavale and A homogeneous twophase Plausibly agreed the experiment, Singhal 2001 approach, with RayleighPlesset with robustness and stability in equations. numerical convergence behaviors Lee et al. 2001a,b Tested the inducer performance, Forwardswept blade was more inlet pressure signals and event resistive to vortex cavitation than characteristics the conventional one Medvitz et al. Homogeneous twophase RANS Obtained qualitative performance 2002 model trends with experimental data in offdesign flow, blade cavitation and breakdown Duttweiler and Surge instability on a cavitation Cyclic oscillation of the attached Brennen 2002 propeller in water tunnel blade cavities Friedrichs and Highspeed camera in a Rotating cavitation was mainly Kosyna 2002 centrifugal pump impeller driven by an interaction of cavity cavitation and the leading edge of the following blades. Hosangadi et al. NS equations with an Well agreed with the experimental 2004 compressible multiphase model data in terms of the head loss and breakdown NSS, and pointed out the performance loss was strongly correlated with cavitation blockage. CoutierDelgosha 3d code FINETURBOTM, with a Obtained good agreement in 2D et al. 2005 timemarching algorithm and low venturi computation, And speed preconditioner for low qualitatively agreed with the Mach number flows experimental measurements of efficiency and cavity visuals in 3D inducer. 3.3 Transport Equationbased Cavitation Model (TEM) Cavitation process is governed by thermodynamics and kinetics of the phase change dynamics in the system. These issues are modeled with the aid of a transport equation with source terms regulating the evaporation and condensation of the phases. In the present study, two different cavitation models, with similar source terms for Eq.(3.3), are presented: Heuristic model (HM). The liquid volume fraction is chosen as the dependent variable in the transport equation (Kunz et al. 2000). The evaporation term mi is a function of pressure whereas the condensation term m' is a function of the volume fraction Cdes va, Min(0, p p) p( (pU2m /2)t, =h. Cproda O (1 _a) (3.7) pit where Cd,, =9.0 x 10 and CPd = 3.0 x104 are empirical constant values, to is the ratio of the characteristic length scale to the reference velocity scale (L /U). Interfacial dynamicsbased cavitation model (IDM). The interfacial dynamics based cavitation model (IDM) was developed by Senocak and Shyy (2004a,b). A hypothetical interface is assumed to lie in the liquidvapor mixture region (Figure 31). By applying the mass and normal momentum conservation equations at the interface, and normalizing with a characteristic time scale t, = L /U, the evaporation mi and condensation tih are given by p laMin(0, p ) P (Vv, Vn)2(p1 )t (3.8) S (1 a)Max(O, p ) (V"" VT"n) 2 (P1 P )t with vapor phase normal velocity, the interface normal, and the interfacial velocity V=un, n=  Va'1 1f pA (3.9) ,n v,n, f =0.9 1 PL wherefis found by computational satisfaction and a value of 0.9 is used. For steady state, we have V,,, = 0.0. This approach will be called original IDM. Vn V Figure 31. Sketch of a cavity in homogeneous flow Regarding the above original IDM, the cavity interfacial velocity was linked to the local fluid velocity in the timedependent computations. Such an approach lacks generality because the interfacial velocity is supposed to be a function of the phase change process. Fundamentally, the interfacial velocity can possibly be estimated more accurately based on the moving boundary computational techniques (Shyy et al. 2004). Here, we estimate the interfacial velocity via an approximate procedure by accounting for the phase transformation process in each computational cell. By integrating Eq. (3.8) through the control volume, we have ihi AV and ir AV. The net interface velocity (the interface velocity relative to the local flow field) becomes Vf"A = AM In n A abs(i AV)abs(i AV) (3.10) '" A A where A is the interface area between vapor and liquid phases. Practically, the control volume face area, Areacy, is projected to the interface normal direction, which can be obtained by taking the gradient of the volume fraction, denoted as S, as shown in Figure 32 for the 2D situation Va  = VctL = fly] n= =n+nJ = Area (3.11) Liquid l,n fLt Interface Areav Vapor Figure 32. Interface vector sketch in a CV Substituting S into Eq. (3.10) leads to vne A abs(i AV) abs(AV) (3.12) V (3S.2) I'" S S The interface velocity includes two parts: flow field local velocity V"LO1 C V= V and net velocity vnt. Then we can have the following derivation ( V 2= [ (et + VLoco ,, (V,) (3.13) Finally, the source terms assume the format, which we call modified IDM SPLLMin(0,p pv) Pv(vt )2(PL Pv)t. (3.14) (1 aL)Max(0, p p ) (Vne )2 (PL PV)tm Recalled that for steadystate condition, we have the relationships of I, = In + VLco" = 0 and v" = _VL, = ,,. In this case, the modified approach is identical to the original IDM. CHAPTER 4 NUMERICAL METHODS The governing equations, presented in the previous chapters, are discretized using a finite volume approach in the present study. In this approach, the flow domain is divided into control volume cells and the governing equations are integrated over each control volume. The main advantage of the finite volume method is that the conservation laws are satisfied locally for each control volume. A nonstaggered grid system is defined at the center of the cell. The decoupling of velocitypressure can be handled by the momentum interpolation method, originally proposed by Rhie and Chow (1983). While, the original momentum interpolation had some problems, such as underrelaxation factor dependent, time step sizedependent, and even checkerboard pressure field, which can be solved using linear interpolation of the two neighbor nodes in the cellface velocity evaluation. A detailed review can be obtained in Yu et al. (2002). In cavitating flow computation, the conventional computational algorithm of singlephase incompressible flow meets severe convergence and stability problems. The situation is improved by using either densitybased method or pressurebased method. Both have been successful to compute turbulent cavitating flows in different configurations with comparable accuracy (Wang et al. 2001, Senocak and Shyy 2004a). However, generally the densitybased method needs preconditioning or artificial density for flows which are largely incompressible (Merkle et al. 1998, Kunz et al. 2000, Ahuja et al. 2001). Hence we choose the pressurebased method (Shyy 1994, Senocak and Shyy 2002a). To take the advantage of noniteration in the timedependent computations, we use the PressureImplicit Splitting Operator (PISO) algorithm other than the Semi Implicit Method for PressureLinked (SIMPLE) method (Patankar 1980). The original PISO method was introduced by Issa (1985) and was modified later to be suitable for large density variation (Thakur and Wright 2002, Thakur et al. 2004). Senocak and Shyy (2002a) further extend it to cavitating flows by addressing the large density jump between phases. 4.1 PressureBased Algorithm The pressurebased algorithm for steadystate computation follows the SIMPLE algorithm (Patankar 1980). The momentum equations can be discretized as AIp = Z Aii Vp(VdP), +b (4.1) where A", A' are the coefficients of the cell center and neighboring nodes from convection and diffusion terms, and Vp, b" are the cell volume and source term separately. Vd is the discrete form of the gradient operator. When there is no source term, the above equation turns into p = H(ib)Dp(VdP)p (4.2) with VP [ /A; 0 Dp= 0 Vp/A' 0 P P The pressurecorrection equation in the pressurebased method has been revised to achieve successful solutions for highly compressible flows (Shyy 1994, Senocak and Shyy 2002a). Generally, the mixture density at the phases interface region has high variations. Here we will illustrate some key computational issues by focusing on the flux terms in the continuity equation, which can be decomposed as pu = (p + p')(i* + i')= p*u* + p*u' + p'u' + p' (4.3) where the starred variables represent the predicted values and primed variables represent the correction terms. And it leads to pp = CP' (4.4) Vd .(pDVdP')+Vd .(Cp*P')= Vd .(p**)+Vd .(CpPDVdP') (4.5) The relative importance of the first and second terms in Eq. (4.3) depends on the local Mach number (Shyy 1994). For low Mach number flows, only the first term prevails; for high Mach number flows, the second term dominates. The fourth term is a nonlinear secondorder term and can be either neglected or included in the source term for stability in early iterations. In the present algorithm, the following relation between density and pressure correction is taken pp =C(1a )pP (4.6) where C is an arbitrary constant and it does not affect the final converged solution. The further details of the model can be obtained from Senocak and Shyy (2002a, 2004a) and Senocak (2002). 4.2 Pressure Implicit Splitting of Operators (PISO) Algorithm for Unsteady Computations In the SIMPLEtype of the pressurebased methods (Patankar 1980, Shyy 1994), the equations are solved successively by employing iterations. In cavitating flow computations, the typical relaxation factors used in the iterative solution process are smaller than the ones used in singlephase flows, and hence smaller time steps are needed to study the cavitation dynamics. Issa (1985) developed the PISO method for the solution of unsteady flows. The splitting of pressure and velocity makes the solution procedure sequential in time domain and enables the accuracy at each time step. It also eliminates the need for severe underrelaxation as in SIMPLE type algorithm. Bressloff (2001) extended the PISO method for highspeed flows by adopting the pressuredensity coupling procedure in allspeed SIMPLE type of methods. Oliveira and Issa (2001) followed the previous PISO work to combine the temperature equation to simulate buoyancydriven flows. Thakur and Wright (2002) and Thakur et al. (2004) developed approaches using curvilinear coordinates with suitability to all speeds. Senocak and Shyy (2002a) further extended this PISO algorithm to enhance the coupling of cavitation and turbulence models and to handle the large density ratio associated with cavitation. The PISO algorithm contains predictor and correction steps. 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 iu = H[Ii], D,(VdP"n1), + u) (4.7) 6t 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 pressure field P* are expected. The discretized momentum equation at this step is written as i* = H[ i], DP(VdP*), + (4.8) P9 Subtracting Eq. (4.7) from Eq. (4.8), we have ui = Dp(VdP'), (4.9) If the density field depends on the pressure field, such as in high Mach number flows or in cavitating flows, the density is corrected using the pressurebased method p = p1 + P, P, = CP' (4.10) The discretized continuity equation written for the new velocity field and density field becomes n 1 P* P" V, +A[p*u** iS ] = 0 (4.11) St Combining Eq. (4.9), (4.10) and (4.11), a first pressurecorrection step equation is obtain below p +A[p" 'D(VdP').S] ],+A[CP'U ] = A[p 'U ] (4.12) 9t To satisfy the mass conservation, the second corrector step is conducted to seek a new velocity field, '*** and pressure field P** iu = H[**] D(Vd**)p + P (4.13) Subtracting Eq. (4.9) from Eq. (4.13) leads to the correction term S= ii + H[i** u ] D (VdP")p (4.14) The corrected density field leads to p =pp + p = p + p' +p (4.15) The second pressurecorrection step can be derived from Eq.(4.14), Eq. (4.15) and the mass continuity equation to reach the below format + p A[p D(VdP").iSCf] + A[CP"U** ] = A[p*U ] 9t (4.16) A[p*H(ti** i*) ). Sf ]p Then, by solving the above predictor and correction steps coupled with the cavitation model, which was formulated in Chapter 3, and turbulence closures given in Chapter 2, the solution procedure for turbulent cavitating flow computations is accomplished. 4.3 SpeedofSound (SoS) Numerical Modeling As mentioned in chapter 1, the harmonic speedofsound in the twophase mixture is significantly attenuated. Therefore the multiphase flow fields are characterized by widely different flow regimes, such as incompressible in pure liquid phase, low Mach compressible in the pure vapor phase, and transonic or supersonic in the mixture. Consequently, an accurate evaluation of speedofsound is necessary and important. From Eq. (4.6) and the definition of speedofsound, the relation between C, and the speedofsound is ap 1 C, = () = (4.17) OP c2 In highspeed flows, the exact form of the speedofsound can be computed easily from the equation of state. However, in cavitating flows, computation of the speedof sound is difficult. Each transport equationbased cavitation model defines a different speedofsound as a result of a more complex functional relationship. In the literature, there have been theoretical studies on defining the speedofsound in multiphase flows (Wallis 1969). Onedimensional assumption and certain limitations are typical in these studies. These definitions do not necessarily represent the actual speedofsound imposed by the cavitation models of interest in this study. They can only be an approximation. On the other hand, the fundamental definition of speedofsound as given in Eq. (4.17) could be useful, and the path to compute the partial derivative is known. From these arguments, it is clear that the computation of the speedofsound in cavitating flows is an open question. Due to the lack of a dependable equation of state for multiphase mixtures modeling sound propagation, Senocak (2002) and Senocak and Shyy (2003) present two numerical forms of speedofsound, called SoS1 and SoS2, and showed that the SoS2 was more likely to produce the correct unsteady behavior in unsteady simulations. In the present study, these two different definitions for the speedofsound are further investigated in the pressuredensity coupling scheme. The SoS1 is the previous pressuredensity coupling scheme with an order of 1 constant coefficient C as ap 1 SoS1: C, = (), = C(1 a) (4.18) The SoS2 is based on an approximation made to the fundamental definition of speedofsound. It is assumed that the path to compute the partial derivative is the OP mean flow direction (4) other than the isentropic direction (s), because the details of thermodynamic properties are not known and the entropy can not be directly computed. This definition is referred to as SoS2 in the rest of the study and given below SoS2: C = ( = Ap ( 1 (4.19) O QP A c2 AP \  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. Since the above approximation in SoS2 is based on the mean flow direction, to generalize it, an averaged form is adopted. It will be referred as SoS2A, taking the format as C coffu P+1,j,k P,1i,j,k + coff Pi,j+1,k P,j1,k C =coffu  +coffrv   P,+1,j,k P,1,],k P,,j+1,k P,,j1,k SoS2A: (4.20) SPl ,j,k+l Pl,j,k1 +coffw   coff ,j,k+ P,j,k 1 with the velocity weight coefficients coffu = Ujk O ,coff = k f ,k + jk + Wk U ,k + ,,,k + W,,,k C &Wlj,k coffw =  ,k + wV,,k U^Jk + V,,k + Wl,,,k CHAPTER 5 ASSESSING TIMEDEPENDENT TURBULENT CAVITATION MODELS 5.1 Cavitating Flow through a HollowJet Valve There are serious implications on the safe and sound operation in flowcontrol valve cavitation phenomena. A limited number of experimental studies have been published on this topic, such as those by Oba et al. (1985) and Tani et al. (1991a,b). In addition, Wang (1999) used highspeed cameras and Laser Light Sheet (LLS) to observe the cavitation behavior in a hollowjet valve under various cavitation conditions and for different valve openings. However, to the best of our knowledge, to date, no comprehensive numerical study has been done in this respect. Furthermore, complex geometries and inaccessible regions of occurrence restrain experimental investigations in cavitation. Hence, we investigate the capability of transport equationbased cavitation models to predict incipient level cavitation. As documented in Wang (1999), Figure 51 shows the geometry and the main configurations of the valve. A key component is the needle, used to control the flow rate by moving to different location in the xdirection. A cylindrical seal supports the needle. There are six struts supporting the cylinder in the pipe center, which are called splitters. The gear is used to control the needle position moving through the xaxis. Figure 52 illustrates the computational domain in selected planes according to the geometry. A multiblock structured curvilinear grid is adopted to facilitate the computation. Figure 52A shows the configuration on the XY plane, and Figure 52B from the YZ plan. Figure 52C shows the plane's location, and Figure 52D shows the boundary conditions in the computations. In the present study, the splitter thickness is neglected and its shape is considered to be rectangular. The Reynolds number is 5x 105 and the cavitation number is 0.9, with the density ratio between the liquid phase and vapor phase p, l/P being 1000 in the water, and the valve opening is 33%. Two investigations are conducted here using the original interfacial dynamics based cavitation model (IDM) (Section 3.3), using the LaunderSpalding model (LSM) (Section 2.2) as turbulence closure. We study steady and unsteady computations. In steadystate, we adopt SoS2. In timedependent simulations, we further examine the two different SoS impact, including SoS1 and SoS2 (Section 4.3), which have been firstly investigated by Senocak and Shyy (2004b) in cavitating flows. For the timedependent computation, the steady singlephase turbulent flow, without considering cavitation, is computed and then the solution is adopted as the initial condition for the unsteady cavitating turbulent flow. The results and discussion are presented in the following. Figure 51. Valve geometry. (1) Splitter. (2) Cylinder. (3) Plunger. (4) Needle. (5) Needle seal overlay. (6) Seal seat inlay. (7) Stroke. (8) Ventilation duct. (9) Gear bZ Upper boundary I :net Splitter outlet Lower boundary alve Pipe center A B Y Lz Noslip BSplitterary c planes tead Noslip s lioutlet Inlet / Noslip Noslip Middle Slis p la n e S lipD C D Figure 52. Computational domains and boundary conditions. A) Computational domain in XY plane. B) Computational domain in YZ plane. C) Planes location. D) Boundary conditions 5.1.1 Steady and Unsteady Turbulent Cavitating Flows Figure 53 shows the density distributions for the steady situation. The cavity is located at the valve tip. The comparison shows that the cavities on the splitter plane and the middle plane are slightly different in size. Figure 54 shows how the cavity shape and location vary with time. For timedependent computations, the cavity fluctuates quasiperiodically. At non dimensional time t* = t/t =0.4, the cavity is the biggest with smallest density, see Figure 54A. As time passes, the cavity size reduces to the smallest in Figure 54(B, C) at time t*=0.6 and t*=0.8, respectively. After one cycle it reaches the maximum again (Figure 5 4D at time t*=1.2s). Both steady and unsteady results exhibit the cavity at the needle tip (Figure 53 and Figure 54). The steady case corresponds to a single instantaneous result in the timedependent solution (Figure 53B and Figure 54C). The timedependent results are qualitatively consistent with the experiment in size and shape (Figure 54). Unfortunately, there is insufficient experimental information to ascertain the time dependent characteristics in detail. The experiment observed that the cavities incept, grow, then detach from the needle tip and transport to the downstream periodically, which is clearly shown in Figure 54E. However, as already discussed in Senocak (2002), with the current combination of turbulence and cavitation models, the detachment of the cavity is not captured, possibly due to the representation of the turbulence via a scalar eddy viscosity. Henceforth, the issue is under further investigation in following by means of adopting different turbulence models. 5 0.99 4 0.98 5 0.99 3 0.96 4 0.98 2 0.9 3 0.96 1 0.8 2 0.9 1 0.8 5 5 4 4 Needle tip Needle tip A B Figure 53. Density contour lines of the steady state solution (The mixture of vapor and liquid inside the outer line forms the cavity), original IDM with SoS1, LSM. A) At splitter plane. B) At middle plane 57 5 099 5 099 4 098 4 098 3 096 3 096 2 09 2 09 1 098 1 08 5 ~ 4 4 Needle Needle A B F5 099 5 099 4 098 4 098 3 096 3 096 2 09 2 09 1 08 1 08 3/ 3 Needle Needle C D Valve tip Figure 54. Middle section density contours at different time instants (The mixture of vapor and liquid inside the outer line forms the cavity), original IDM with SoS2, LSM. A) Time t*=0.4. B) Time t*=0.6. C) Time t*=0.8. D) Time t*=1.2. E) Experimental observations from Wang (1999): cavity at needle tip and cavitation aspects around the needle To further demonstrate the cavityinduced quasiperiodic characteristics of the flow field, Figure 55 highlights the time evolutions at selected locations. Figure 55A shows the locations of the points selected. The pressure is at the middle points, and the density is near the bottom boundary, just one point away from it. From Figure 55B, at the places away from the cavity, A, B and D, the density is constant since only the liquid phase exists there, and is quasiperiodic inside the cavity at Point C. On the other hand pressure oscillates in the whole domain except on the inlet plane, where the flow condition is fixed (Figure 55C). C D A B, A 095 7  SLocation A 6  09 Location B 5 Location A SLocation C 4 Location B G Location D V Location C 085 C o 3 Location D 2  08 2 08 0 75 07 2 0 1 2 3 4 0 1 2 3 4 SB t C Figure 55. Time evolutions at different locations, original IDM with SoS2, LSM. A) Samples locations. B) Density. C) Pressure Figure 56 presents the flow structure on the middle plane. For the steady case, the flow fields of both single phase (without invoking the cavitation model in the course of computation) (Figure 56A) and cavitating flows are almost identical (Figure 56B). It 59 indicates that in the present case, the detailed cavitation dynamics does not exhibit substantial influence on the overall flow pattern. For the timedependent case, at different time instants, t*=0.4 and t*=0.8 (Figure 56C), the flow field around the needle remains largely the same. The comparison between the steady and unsteady case indicates that there is not very much difference in the flow pattern. As expected, there are two recirculating zones: one behind the needle and the other one downstream around the splitter region, which is located at about x=3.44.2. Compared with the experimental observation in Figure 56D (the cavitating flow structures behind the needle, and in the splitter region), the present results are in general agreement. 1 0  0 0 3 35 4 45 5 3 35 4 45 5 55 0 0  3 35 4 45 5 3 35 4 45 5 X X C Figure 56. Projected 2D streamlines at middle plane and experimental observation. A) Steady singlephase flow, LSM. B) cavitating flow, original IDM with SoS2, LSM (Right). B) Unsteady cavitating flow at time t*=0.4 (Left) and t*=0.8 (Right). D) Flow pattern from Wang (1999), behind valve about x=3.03.4 (Left) and at splitter region about x=3.44.2 (Right) Flow  Figure 56. Continued 5.1.2 Impact of Speedofsound Modeling We further examine the behavior of the different handling of speedofsound (Chapter 4), that is, SoS1 and SoS2. Figure 57 shows the variation in cavity shapes and location for LSM at maximum and minimum fraction instants. Compare the half period between maximum and minimum instant, the cavity oscillation for SoS1 is distinguished by a significantly higher frequency and moderately higher amplitude than for SoS2 (Figure 57A, B). These observations, as also discussed later, are a precursor to the significant impact produced by the speedofsound definition on the flow time scales. The qualitative shape and location of the cavity for both SoS1 and SoS2, however, are similar to the experimental observations (Figure 56D). Senocak and Shyy (2003,2004b) pointed out that SoS2 performed better than SoS1 in simulating the oscillatory behavior of the cavitating flow in the convergentdivergent nozzles, because SoS2 successfully produced a quasisteady solution at high cavitation number and a periodic solution at low cavitation number. In Figure 58, SoS2 and SoS1 behave similarly in this computed 61 case, and they both are nearly constant when away from the pure phase. These qualitatively agree with the theoretical analysis (Figure 12). Also, this helps to explain why the cavity shapes of the two different SoS models are qualitatively close. One thing needs to be pointed out here that the volume fraction is very close to the pure liquid phase. Hence the linear approximation SoS1 can have similar behavior as SoS2. 5 099 4 098 3 096 2 09 1 08 5 099 4 098 3 096 2 09 1 08 5 4 5 099 4 098 3 096 2 09 1 08 // ^ 5 099 1 4 098 S3 096 2 09 1 1 08 ^^ ^4 6 Figure 57. Middle section density contours of different SoS at different instants (the mixture of vapor and liquid inside the outer line forms the cavity), original IDM, LSM. A) SoS2, maximum fraction at t*=0.32 (Left) and minimum fraction at t*= 0.77 (Right). B) SoS1, maximum fraction at t*=5.92 (Left) and minimum fraction at t*=5.86 (Right) Lack of substantial experimental data prevents us from making more direct comparisons at this stage. Various experimental studies have shown that, besides oscillating, the cavities incept, grow, detach from the needle tip and get periodically 3// transported downstream (Figure 54E). The periodic detachment and inception of cavities is, however, difficult to capture through currently known turbulence and SoS models (Wang et al. 2001, Senocak and Shyy 2003, CoutierDelgosha et al. 2003, Wu et al. 2003 a). Consequently, fundamental understanding on the above modeling aspects, which is discussed in the following sections, is imperative to encounter the above limitation. The following discussions serve to probe the sensitivity of the various modeling concepts, guided by our qualitative, but incomplete, insight into the fluid physics. The implications of speedofsound (SoS) definition are quantified by a correlation study (Table 51 and Table 52) and spectral analysis (Figure 510) on a series of pressuredensity time history at various nodes in cavitation vicinity. The Pearson's correlation (r) between pressure and density is calculated as N Y[(X, )(Y )] r= I ; X =P, Y= p (5.1) V x, (XrY)2 Table 51 and Table 52 indicate that within similar variation limits of liquid volume fraction, SoS1 has a broader range of pressuredensity correlation coefficients than SoS2. Furthermore, SoS2 consistently exhibits much stronger pressuredensity coefficients than SoS1. The distinct effects produced by the speedofsound definition are further corroborated by the dramatic time scale differences observed from the time history and FFT plots of pressure and density (Figure 59 and Figure 510A, B, C). SoS1 plots clearly indicate dominance of a single highfrequency compressibility effect unlike SoS2, which are characterized by a wider bandwidth for pressure and density, and a lower dominant frequency. It is worth noting here that although the plots in Figure 59(B, 63 C) are plotted against a smaller time range, we have used a convincingly long time history for our analysis. 70 60 0 40 0 U, 0 09 095 liquid fraction Liquid fraction Figure 58. Two different SoS in hollowjet valve flow, original IDM, LSM Table 51. Timeaveraged liquid volume fraction v/s pressuredensity correlation at multiple points inside the cavity, original IDM with SoS1, LSM Time averaged liquid volume Pearson's correlation (r) between fraction (al) pressure and density (SoS1) 0.862 0.674 0.882 0.09 0.936 0.093 0.970 0.045 0.978 0.077 0.983 0.096 0.986 0.103 0.988 0.12 Mean(r)= 0.13; SD(r)= 0.9 Table 52. Timeaveraged liquid volume fraction v/s pressuredensity correlation at multiple points inside the cavity, original IDM with SoS2, LSM Time averaged liquid volume Pearson's correlation (r) between fraction (al) pressure and density (SoS2) 0.858 0.733 0.868 0.106 0.934 0.106 0.968 0.46 0.977 0.47 0.982 0.412 0.985 0.412 0.987 0.412 Mean(r)= 0.42; SD(r) = 0.1 SLocation A Location B Location C SLocation D . E B Figure 59. Pressure time evolutions of different locations. B)SoS1. C) SoS2 Location A Location B Location C Location D 2 0 1 2 tS SoS, original IDM, LSM. A) Samples I I I IA Figure 510. Time evolution and spectrum of pressure and density of two SoS definitions at a point at the cavitation vicinity, original IDM, LSM. A) Pressure (Left) and density (Right) time history. B) Spectral analysis on pressure (Left) and density (Right, SoS1. C) Spectral analysis on pressure (Left) and density (Right), SoS2 s " B I I I3 2 25I C Figure 510. Continued 5.2 Turbulent Cavitating Flow through a ConvergentDivergent Nozzle To study the performance of the filterbased model (FBM, Section 2.3), and the modified IDM (Section 3.3) in timedependent computations, the turbulent cavitating flow through a convergentdivergent nozzle is investigated. Stutz and Reboud (1997, 2000) have studied the unsteady cavitation formed in a convergentdivergent nozzle that has a convergent part angle of 180 and a divergent part angle of 80. The experimental Reynolds number is Re = 106 3 3 x 106 based on the reference velocity and cavity length. Cavitation formed in this nozzle is described as "unsteady and vapor cloud shedding". Senocak and Shyy (2003, 2004b) conducted the numerical simulations using IDM coupled with the LSM. They found that under cavitation number c = 1.98, the cavity length L,, matched the experiments (Stutz and I $a 0 I 2D 5 rhapy  I I I I I ILm q. mm U l I" ir m d m, J .. b Reboud 1997, 2000), which did not specify the cavitation number. Therefore, to be consistent, we use c = 1.98 in the present study. First, we apply the original IDM with SoS2 for the simulations to assess the FBM performance. To test the grid and filter size dependency, two grids are adopted including a coarse grid with A = 0.5Lo and a fine grid with A = 0.25L av. The timeaveraged eddy viscosity contours of two grids are given in Figure 511. The FBM gives effective viscosities for both grids about half an order of magnitude lower than that produced by the LSM. The differences of the timeaveraged eddy viscosity are noticeable between the two grids. The largest viscosity on the fine grid is lower and its location shifts to the upstream, compared with Figure 511 (A, B) left and right columns correspondingly. However, the time averaged vapor volume fraction does not show significant difference (Figure 512A, B). The timeaveraged uvelocity profiles have the similar trend. Hereafter, all the results are based on the fine grid with A = 0.25Lv Figure 513 shows the timeaveraged velocity and vapor volume fraction profiles within the cavity at four different sections, using two different turbulence models. The boundaries of the cavitating region, from experiments, LSM and FBM, are also included. Although the numerical results of the cavitating boundaries are about 5% higher than that from the experiments, the computations capture the main cavity body and the overall trends are agreeable (Figure 513B). The velocity profiles qualitatively agree with the experimental data, especially in the core of the reverse flow (Figure 513A). Comparing the LSM and FBM, both the timeaveraged velocity and volume fraction show marginal difference. However, the instantaneous velocity and volume fraction profiles display significant differences (Figure 514). The current results show an autooscillating cavity. Here only the maximal and minimalvolume fraction time instants are presented. The maximal volume fraction of the FBM is bigger than that of the LSM. On the contrary, the minimal volume fraction is much smaller by comparing Figure 514A and Figure 514B. The reverse velocity profiles show the same trend. The increase in fluctuations is due to the FBM reducing the eddy viscosity (Figure 511). The observation reveals that the time averaging process can be misleading and not suitable as the only indicator for performance evaluation, which is similar to the previous study by Wu et al. (2003b). I 4111E 11 1005 U U 604 E04 005 005 0 41 0 4904 15 LSM o i,. FBM 05 0 5 1 15 05 0 05 1 15 x (m) x(m) S25617E04 7 2 4 56E04  4 3541E 044 1E04 2 4191E04191E LSM FBM 0 15 0 15  05 0 D5 1 15 0 5 o 05 1 15 x(m) x(m) B Figure 511. Timeaveraged eddy viscosity contours of different grids, original IDM with SoS2, = 1.98. A) Coarse grid with A= 0.50LZav. B) Fine grid with A= 0.25L A = 0.25L, Vapor volume fraction, coarse grid Vapor volume fraction, fine grid Vapor volume fraction, Exp. data Boundary of cavitating flow, coarse grid Boundary of cavitating flow, fine grid 25% 25% 25% 25(%) 25% Vapor volume fraction, coarse grid, A=0.5L_ Vapor volume fraction, fine grid,A=0.25L_ 15 Vapor volume fraction, Exp. data Boundary of cavitating flow, coarse grid Boundary of cavitating flow, fine grid 25% 25% 25% 10 25% 5 n 0 20 40 60 80 0 20 40 60 80 X (mm) X (mm) A B Figure 512. Timeaveraged vapor volume fraction comparisons of different grids, original IDM with SoS2, a = 1.98. The vertical scale is the distance from the wall. Experimental data are from Stutz and Reboud (1997&2000). A) LSM. B) FBM Velocity, Standard Vapor volume fraction, Standard Velocity, filterbased Vapor volume fraction, filterbased 15 A Velocity, Exp. data 15 A Vapor volume faction, Exp. data Boundary of cavitating flow, Exp. Boundary of cavitating flow, Exp. Boundary of cavitating flow, Standard Boundary of avtating flow, Standard Boundary of cavitating flow, Filterbased U (m/s) Boundary of cavitating flow, Filterbased (%) 10 10 25% 25% 10 25% 10 A 10 A 25% 5 E 7 A E X (mm) A X (m) A B Figure 513. Timeaveraged comparisons of different turbulence models on fine grid with A = 0.25L original IDM with SoS2, a = 1.98. The vertical scale is the distance from the wall. Experimental data are from Stutz and Reboud (1997&2000). A) Uvelocity. B) Vapor volume fraction An illustration that the FBM can induce noticeably stronger flow oscillations is presented in Figure 515. It is observed that timeaveraged pressure contours and streamlines are similar at the cavity zone, although a recirculation region produced by the FBM exists at downstream near the lower wall (Figure 515A). However, the I // ii 69 instantaneous pressure contours and velocity field produced by the two turbulence models are distinctly different. The FBM yields wavy flow structures, which induce the autooscillation of the flow fields in Figure 515(B, C). The experimental measurements (Figure 514) are confined in the throat region, which have missed the large scale unsteadiness of the flow field.  Velocity, minimal instant Velocity, maximal instant 15 Boundary of cavitating flow, Minimal Boundary of cavitating flow, Maximal 10 10 L /j 'I 10 S 10 20 30 40 50 X(mm) 60 70 80  Velocity, minimal instant Velocity, maximal instant 5 Boundary of cavitating flow, Minimal S Boundary of cavitating flow, Maximal 10 / 10 * I II u (m/s) 10 //  Vapor volume fraction, minimal instant Vapor volume fraction, maximal instant Boundary of cavitating flow, Minimal Boundary of cavitating flow, Maximal (y) 25% 25%  Vapor volume fraction, minimal instant Vapor volume fraction, maximal instant  Boundary of cavitating flow, Minimal Boundary of cavitating flow, Maximal ^ (o) 25% 0 10 20 30 40 50 60 70 80 _0 10 20 30 40 50 60 70 80 X (mm) X (mm) Figure 514. Instantaneous profiles on fine grid with A = 0.25Lc, original IDM with SoS2, a = 1.98. The vertical scale is the distance from the wall. A) LSM: u velocity and vapor volume fraction. B) FBM: uvelocity and vapor volume fraction I I P j Pressure contours: color S0.842 0.1  LSM 0 05 x (m) Pressure contours: color 0.842 LSM 05 x(m) Pressure contours: color 0.842 Pressure contours: color *0.842 LSM 11.26 05 S(m) Figure 515. Pressure contours and streamlines comparison of two turbulence models on fine grid, original IDM with SoS2, cr = 1.98. A) Timeaveraged. B) At maximalvolume fraction instant. C) At minimalvolume fraction instant Obviously, both turbulence and cavitation models directly affect the cavity dynamics such as the shedding pattern and frequency. Hence in this section, we further investigate the implication of the cavitation model in the context of the interfacial velocity estimated, as discussed in Section 3.3 which is called the modified interfacial dynamicsbased cavitation model. The speedofsound (SoS) is based on the SoS2A, formulated in Section 4.3. The turbulence closure is employed the original Launder Spalding model (LSM) described in Section 2.2. The detailed results are presented below. Figure 516 gives the time history of the pressure at a reference point near the inlet. The pressure of the original IDM, which adopts an empirical factor to assign the interfacial velocity based on the local fluid velocity in Section 3.3 Eq. (3.9), exhibits much smaller variation than that computed with the modified IDM. As given in table 53, although both IDMs underpredict the Strouhal number, the dynamically adjusted interfacial velocity embodied by the modified IDM yields better agreement compared with the experimental measurement. Selected time evolutions of the cavity shape and the associated flow structures are presented in Figure 517 with the available experimental observation. The original IDM demonstrates an attached periodic cavity, while the modified IDM produces cavity break up and detachment, in manners qualitatively consistent with the experimental observations. The cavity breakup and collapse of the modified IDM also helps to explain the Cp behavior shown Figure 516. Furthermore, because of the shedding of the cavity, the modified IDM better predicts the vapor volume fraction profiles than the original IDM, with less difference in the velocity profiles between these two cavitation models (Figure 518). Table 53. Comparisons of Strouhal number of original and modified IDM with SoS2A Original IDM, SoS2A Modified IDM, SoS2A Experiment 0.07 0.13 0.30  Modified IDM Original IDM HP y 1 2 3 4 5 Figure 516. Pressure evolutions of original reference point and modified IDM with SoS2A at a Figure 517. Cavity shape and recirculation zone during cycling of original and modified IDM. The experimental observation is adopted from Stutz and Rebound (1997, 2000) .1,i  Figure 517. Continued X (mm) X (mm) A B Figure 518. Timeaveraged volume fraction and velocity comparisons of original and modified IDM. The vertical scale is the distance from the wall. Experimental data are from Stutz and Reboud (1997, 2000). A) Vapor volume faction. B) u velocity 5.3 Turbulent Cavitating Flow over a ClarkY Hydrofoil The capability of the modified IDM (Section 3.3) and turbulence models are further investigated in unsteady cavitating flows over a ClarkY hydrofoil, assessed by experimental data from Wang et al. (2001). To facilitate the performance evaluation, the generalized treatments of the SoS2, including SoS2, SoS2A described in Section 4.3, are coupled. For the computational setup (Figure 519A), the computational domain and boundary conditions are given according to the experimental setup. The ClarkY hydrofoil is located at the water tunnel center. Two anglesofattack (AoA) considered are 5 and 8. The hydrofoil chord length is c and the hydrofoil leading edge is 3c away from the inlet. The two important parameters are the Reynolds number and the cavitation number, which is based on inlet pressure P. and vapor pressure P1 with inlet velocity U. Re Uc 7 x 105 (5.2) V P P S= (5.3) pU /2 The filter size in the present study is chosen to be larger than the largest grid scale employed in the computation, and is set to be A = 0.08c. Computations have been done for two AoA and several cavitation numbers. Specifically, for AoA=5 under four flow regimes: nocavitation (o = 2.02), inception (o = 1.12), sheet cavitation (o = 0.92), and cloud cavitation (o = 0.55), and AoA=8 under three cavitation numbers: nocavitation (o = 2.50), sheet cavitation (o = 1.40), cloud cavitation (o= 0.80). All cases above are at the same Reynolds number Re =7x105. Two different turbulence models, LSM and FBM, have been employed to help probe the characteristics of cavitation and turbulence modeling interactions, coupled with different SoS numerical models. The outlet pressure is fixed and adjustment of the vapor pressure is needed to be consistent with the prescribed cavitation number in each case. Grid sensitivity analysis. To investigate the grid dependency, two grids are adopted in the computation: coarse grid and fine grid. The grid blocks and numbers of the coarse grid are shown in Figure 519B. The fine grid has 60% more number of nodes than that of the coarse grid in the vertical direction while maintaining the same distribution in the horizontal direction. Two different cavitation numbers, o = 2.02 and a = 0.55, both at AoA=5 have been investigated. Overall, the solutions on both grids are in good agreement. The timeaveraged uvelocity and vvelocity profiles with o = 2.02 and a = 0.55 using the LSM are shown in Figure 520. The results based on the FBM are of similar nature and will not be repeated. Hereafter, to reduce the cost of timedependent computations, we use the coarse grid in the computations. Visualization of cavity and flow field. First, we focus on the LSM results to analyze the performance of the IDM. Figure 521 and Figure 522 show the time averaged flow structure and cavity shape under varied cavitation numbers. With no cavitation, the flow field is attached without separation for both AoA=5 and AoA=8 (Figure 521A left column and Figure 522A left column). This is consistent with the experimental observation. When cavitation appears, the density will change by a factor of 1000 between liquid and vapor phases. Consequently, there is a drastic reduction in the amount of mass inside the cavity, and a contraction of the fluid flow behind the cavity. With the reduction of the cavitation number, the cavity and recirculation zone become bigger. At cloud cavitation regime, the cavity experiences shedding, causing multiple recirculating flows (Figure 521D left column and Figure 522C left column). Compared to the experimental data for both AoAs (Table 54 and Table 55), the predicted cavity sizes demonstrate qualitatively consistent rending, albeit generally overpredicted. For incipient cavitation, the experiment observed the recurring formation of hairpin type cavitating vortex structures, which are not attached to the solid surface (Wang et al. 2001). This type of flow structure is not captured in the computation. The timeaveraged flow structures associated with sheet and cloud cavitation, on the other hand, seem to be reasonably captured computationally. Furthermore, the timeaveraged outcome of employing both LSM and FBM seems compatible. NoSlip e i Hydrofoil: NoSlip _E 0 10c NoSlip A Coarse Grid Block 2:70x50 Block 4:70x50 Block 1: 80x50 Block 1: 70x50 /Block 3: 70x50 Block 5: 80x50 B Figure 519. ClarkY geometry sketch and Grid blocks. A) Geometry configuration and boundary conditions, c is the hydrofoil chord. B) Grid blocks and coarse grid numbers 77 A Experiment Fine grid O Coarse grid 04 I o 02  4 1.0 1. 0_ulu  02 04 06 08 1 12 X/C A Experiment  Fine grid Coarse grid 04 02 S1.. 01.. 1.0 002 04 06 08 1 XIl A A Experiment Fine, LSM 4 Coarse, LSM 04  A A A A Aa S 1.01 u/U 1.0 0< Uv/U 1.0 02 04 06 08 1 12 X/c A Experiment a Fine, LSM 04 Coarse, LSM A A 02 I *A A , S.VU A 0 02 04 06 08 1 XC B Figure 520. Grid sensitivity of the timeaveraged u and vvelocity, LSM, AoA=5. Experimental data are from Wang et al. (2001). A) Nocavitation. B) Cloud cavitation c = 0.55 Table 54. Timeaveraged cavity leading and trailing positions of different turbulence models (Taking aL = 0.95 as cavity boundary), modified IDM with SoS2A, AOA=5 degrees Position Inception = 1.12 Sheet cavitation a = 0.92 Cloud cavitation a = 0.55 LSM FBM Exp. LSM FBM Exp. LSM FBM Exp. Leading 0.13 0.13 0.13 0.13 0.23 0.16 0.17 0.22 Trailing 0.62 0.57 0.64 0.64 0.68 1.01 1.04 0.90 Table 55. Timeaveraged cavity leading and trailing positions of different turbulence models (Taking aL = 0.95 as cavity boundary), modified IDM with SoS2A, AoA=8 degrees. Position Sheet cavitation a = 1.40 Cloud cavitation a = 0.80 LSM FBM Exp. LSM FBM Exp. Leading 0.094 0.098 0.15 0.12 0.13 0.15 Trailing 0.53 0.52 0.56 1.07 1.10 0.84 Figure 521. Timeaveraged volume fraction contours and streamlines of different turbulence models, AoA=5. A) Nocavitation. B) Incipient cavitation. C) Sheet Cavitation c = 0.92. D) Cloud cavitation o = 0.55 B C Figure 522. Timeaveraged volume fraction contours and streamlines of different turbulence models, AoA=8. A) Nocavitation. B) Sheet cavitation a = 1.40. C) Cloud cavitation a = 0.80 The temporal evolution of the computed and experimentally observed flow structures with cloud cavitation, under two AoAs, are shown in Figure 523 and Figure 5 24. Figure 523A left column shows the time sequences of flow structures predicted by LSM at AoA=5. The corresponding flow predicted by the FBM is shown in the right column. The experimental visual image is shown in Figure 523B. The flow structures at AoA=8 are shown in Figure 524. Both the computations and the experiment indicate that as the AoA increases, the cavity exhibits a more pronounced recurrence of the size variation. The FBM predicts stronger timedependency than the LSM. We will further investigate this aspect in the following discussion. To help elucidate the main features of the cavity dynamics from both numerical and experimental studies, we show in Figure 5 25 three stages of the cavity sizes, for cloud cavitation c = 0.55 at AoA=80. The figure demonstrates that the numerical simulation is capable of capturing the initiation of the cavity, growth toward trailing edge, and subsequent shedding, in accordance with the qualitative features observed experimentally. 80 t=0.00 t=3.36 ms t= 10.1ms t=11.8 ms t=15.1 ms t=20.2 ms A t=0.0 t=2.2 ms t=4.4 ms t=6.7 ms t=8.9 ms t=ll.1 ms B Figure 523. Time evolutions of cloud cavitation c = 0.55, AoA=5. A) Numerical results of different turbulence models. B) Side views of the experimental visuals from Wang et al. (2001) Figure 524. Time evolutions of cloud cavitation c = 0.80, AoA=8. A) Numerical results of different turbulence models. B) Side views of the experiment visuals from Wang et al. (2001) LSM. t=24.5 t=22.2 ms Figure 524. Continued t=)U ms Flow Figure 525 Cavity stage comparisons, cloud cavitation c = 0.80, AoA=8: FBM (Left) and Experiment (Right)0. Experimental visuals are from Wang et al. (2001). A) Early stage: cavity formation. B) Second stage: cavity growth toward trailing edge. C) Third stage: cavity breakup and shedding Velocity profiles and lift/drag coefficients. The mean horizontal u velocity and vertical vvelocity of the flow field are illustrated in Figure 526 (AoA=5) to Figure 527 (AoA=8). The timeaveraged velocity profiles are documented at 6 chordwise locations, 0%, 20%, 40%, 60%, 80%, and 100% of the leading edge, under different cavitation numbers. With no cavitation, the numerical results agree well with the experiment, and the results of the two turbulence models are virtually identical (Figure 526A and Figure 527A). With the cavitation number decreasing, the differences between prediction and measurement become more substantial, especially at the cavity closure region. In fact, recirculation breaks up the cavity and tilts the rear portion of the cavity more upward, which makes the cavity thicker in comparison to the experiment. With the reduction in cavitation number, it becomes more distinguished in Figure 523 and Figure 524. The cavity tiltup results in overshooting the uvelocity at the closure region (0.8c and 1.0c sections) (Figure 526 and Figure 527). Overall, considering the difficulties in experimental measurement (Wang et al. 2001) the agreement is reasonable. A Experiment  Filterbased LaunderSpalding A Experiment  Filterbased  LaunderSpalding A Expenment  Filterbased LaunderSpalding  FBM, cavity boundary  LSM, cavty boundary A Experiment Filterbased LaunderSpalding FBM, cavity boundary LSM, cavity boundary "xC B Figure 526. Timeaveraged u and vvelocities of two turbulence models, AoA=5. Experimental data are from Wang et al. (2001). A) nocavitation. B) Inception a = 1.12. C) Sheet cavitation a = 0.92. D) Cloud cavitation c = 0.55 04 02 02 12 A Expenment Filterbased SLaunderSpalding ... FBM, cavity boundary LSM, cavity boundary A Experiment  Filterbased LaunderSpalding FBM, cavity boundary LSM, cavity boundary Figure 526. Continued 04 02 04 02 12 