Simulations of Cavitating Flows Using Hybrid Unstructured

(1)

Vineet Ahuja

Research Scientist

Ashvin Hosangadi

Principal Scientist

Srinivasan Arunajatesan

Research Scientist

Combustion Research and Flow Technology, Inc. (CRAFT Tech), P.O. Box 1150, Dublin, PA 18917

Simulations of Cavitating Flows Using Hybrid Unstructured

Meshes

A new multi-phase model for low speed gas/liquid mixtures is presented; it does not require ad-hoc closure models for the variation of mixture density with pressure and yields thermodynamically correct acoustic propagation for multi-phase mixtures. The solution procedure has an interface-capturing scheme that incorporates an additional scalar transport equation for the gas void fraction. Cavitation is modeled via a finite rate source term that initiates phase change when liquid pressure drops below its saturation value. The numerical procedure has been implemented within a multi-element unstruc- tured framework CRUNCH that permits the grid to be locally refined in the interface region. The solution technique incorporates a parallel, domain decomposition strategy for efficient 3D computations. Detailed results are presented for sheet cavitation over a cylindrical head form and a NACA 66 hydrofoil. 关DOI: 10.1115/1.1362671兴

Introduction

Numerical simulations of cavitating flows are very challenging since localized, large variations of density are present at the gas/

liquid interface while the remainder of the flow is generally incompressible. Furthermore, the cavitation zone can detach due to the influence of a re-entrant jet and convect downstream. This detachment/collapse of the cavity in the pressure recovery region can lead to strong acoustic disturbances. These disturbances can manifest themselves as high frequency noise that has the potential to considerably alter the acoustic signal/signature and degrade performance. In general, current cavitation models, which have originated from incompressible formulations, deal with these issues by defining ad-hoc closure models relating density and pressure which are not very general. In our paper, we describe a compressible, cavitation formulation which yields the correct acoustical and thermodynamic behavior for multi-phase systems.

The importance of being thermodynamically consistent is that the same cavitation model can be used for simulations from marine propellers to cryogenic pumps. In particular, the system does not require user defined closure models for capturing the gas/liquid interface and is designed to perform efficiently in the nearly incompressible liquid regime. We give a brief review of the literature below and highlight the differences with our formulation.

Cavitation models in the literature may broadly be classified into two categories; interface fitting and interface capturing procedures. Interface fitting procedures explicitly track and fit a dis- tinct gas/liquid interface which is an internal boundary. While this procedure gets around the numerical difficulties of integrating through the interface, its applications are limited to simpler problems where the cavity can be described as a well-defined closed volume of pure gas. Interface capturing schemes, where the gas/

liquid interface is obtained as part of the solution procedure, are more general in their applications and may be applied to both sheet cavitation as well as bubbly cavitation. Here, the thermodynamic issues of integrating through gas/liquid mixtures with density and acoustic speed variations have to be dealt with in order to provide closure to the equation system.

Closure models reported in the literature are generally tailored for specific problem classes and do not address the multi-phase physics in a fundamental fashion. For instance, bubbly cavitation

has been tackled by providing pseudo-density relations derived from various modified forms of the Rayleigh-Plesset equation 共Kubota et al. 关1兴 and Chen and Heister 关2兴兲. Here, no additional equation for the convection of gas species is specified and the formulations make assumptions about the bubble number density in the fluid. On the other hand, for sheet cavitation an ad-hoc pseudo-density relationship based on the local pressure is provided to close the system共Delaunay and Kueny 关3兴 and Janssens et al. 关4兴兲 which can be quite restrictive since it is decoupled completely from the local phase composition.

More generalized cavitation formulations have been presented in recent papers by Merkle et al.关5兴 and Kunz et al. 关6兴. Here, an additional equation for the gas void fraction is solved for and the local mixture density is obtained from the local phase composition. The equations are cast in a compressible-like time marching format with preconditioning to overcome numerical stiffness problems. However, in both cases, the authors implicitly specified a pressure-density relation which leads to an erroneous acoustic speed in the mixture, particularly in the interface region. The non- physical acoustic speeds are clearly not appropriate for unsteady simulations. However, even for steady calculations, this may pose numerical difficulties at the gas/liquid interface since the physical acoustic speed behaves in a nonlinear fashion共as we shall discuss later兲 and has a very small magnitude; the local Mach number at the interface can be large leading to very different numerical characteristics.

The formulation presented in this paper is an acoustically accurate form of the compressible multi-phase equations共Ahuja et al.

关7兴兲 and is an extension of our earlier work in high pressure gas/

liquid systems 共Hosangadi et al. 关8兴兲. The numerical algorithm follows a similar time-marching philosophy 共as that in Merkle et al.关5兴 and Kunz et al. 关6兴兲. However, here the local speed of sound in the two-phase mixture is a function of the local void fraction and mimics the two-phase acoustic speed relationship from classical analytic theory. We note that the system presented here is closed and does not require additional equations to resolve the gas/liquid interface. The model is general and may be applied to both sheet and bubbly cavitation by providing appropriate source terms for gas generation/reabsorption in cavitating regions.

The numerical formulation has been implemented on a hybrid unstructured framework which permits tetrahedral/prismatic cells.

An unstructured framework is particularly suitable for geometri- cally complex systems like marine propellers where the blades are

Contributed by the Fluids Engineering Division for publication in the JOURNAL OFFLUIDSENGINEERING. Manuscript received by the Fluids Engineering Division April 5, 2000; revised manuscript received January 29, 2001. Associate Editor:

Y. Matsumoto.

(2)

skewed strongly. In addition, the ability for local grid refinement in a complex flow field provided a strong motivation for using unstructured grids here. For instance, in cavitation simulations, the region near the gas/liquid interface exhibits strong gradients in flow properties and requires high local grid resolution which can be achieved most economically with a grid adaption procedure.

For efficient computations of large 3D problems, a parallel framework for distributed memory systems has been implemented. De- tails of the numerical formulation are provided in the following section followed by details of the unstructured framework. In the Results section, we discuss details of our simulation for a cylindrical headform and a NACA 66 hydrofoil. The details of the flowfield in the closing region of the cavity and the turbulence characteristics of the wake are examined in depth. Surface pressure comparisons with experimental data are shown for a range of cavitation numbers to validate the numerics.

Multi-Phase Equation System

The multiphase equation system is written in vector form as:

⳵Q

⳵t ⫹

⳵E

⳵x ⫹

⳵F

⳵y ⫹

⳵G

⳵z ⫽S⫹Dv (1)

Here Q is the vector of dependent variables, E, F, and G are the flux vectors, S the source terms and D_v represents the viscous fluxes. The viscous fluxes are given by the standard full compressible form of Navier Stokes equations共see Hosangadi et al. 关9兴 for details兲. The vectors Q, E, and S are given below with a detailed discussion on the details of the cavitation source terms to follow later:

Q⫽

冉

^␳^␳^␳^␳^␳^␳^␳^g^m^m^m^m^m^␾^m^w^u^v^k^␧^g

冊

^E^⫽

冉

^␳^␳^m^␳^␳^␳^␳^␳^gû^m^m^m^m^␾^m²ûwûv^␧u^kuû^⫹P^gû

冊

^S^⫽

冉

^S^S^S⁰⁰⁰⁰^g^k^␧

冊

⁽²⁾

Here,␳mis the mixture density, and␾gis the volume fraction or porosity for the gas phase. Note that in Eq.共2兲 the energy equation is dropped since each phase is assumed to be nearly incompressible thereby decoupling the pressure work term. An additional scalar equation for mixture enthalpy may be solved coupled to this equation set if the temperature effects become important. The mixture density and gas porosity are related by the following relations locally in a given cell volume:

␳m⫽␳g␾g⫹␳L␾L (3)

1⫽␾g⫹␾L (4)

where ␳g,␳L are the physical material densities of the gas and liquid phase respectively.

To modify the system in Eq.共1兲 to a well-conditioned form in the incompressible regime requires a two-step process; an acoustically accurate two-phase form of Eq. 共1兲 is first derived, followed by a second step of time-scaling or preconditioning to obtain a well-conditioned system. We begin by defining the acoustic form of density differential for the individual gas and liquid phase as follows:

d␳g⫽ 1

c_g²d P, d␳l⫽1

c_L²d P (5)

Here c_g is the isothermal speed of sound (⳵P/⳵␳g)_Tin the pure gas phase, and c_Lis the corresponding isothermal speed of sound in the liquid phase, which is a finite-value. The differential form of the mixture density␳m is obtained by differentiating Eq.共3兲, and using the relationship given in Eq.共5兲 to obtain,

d␳m⫽共␳g⫺␳L兲d␾g⫹ 1 c_␾²d P

冉

^c¹␾2⫽␾g

c_g²⫹␾L

c_L²

冊

⁽⁶⁾

Here, c_␾ is a variable defined for convenience and is not the acoustic speed, c_m, in the mixture which will be defined later.

Using Eq.共6兲, Eq. 共1兲 may be rewritten as:

⌫⳵Qv

⳵t ⫹

⳵E

⳵x ⫹

⳵F

⳵y ⫹

⳵G

⳵z ⫽S⫹Dv (7)

where

⌫⫽

¨

1

c_␾² 0 0 0 共␳g⫺␳L兲 0 0 u

c_␾² ␳m 0 0 共␳g⫺␳L兲u 0 0 v

c_␾² 0 ␳m 0 共␳g⫺␳L兲v 0 0 w

c_␾² 0 0 ␳m 共␳g⫺␳L兲w 0 0

␾g

c_g² 0 0 0 ␳g 0 0

k

c_␾² 0 0 0 共␳g⫺␳L兲k ␳m 0

␧

c_␾² 0 0 0 共␳g⫺␳L兲␧ 0 ␳m

©

⁽⁸⁾

and,

Q_v⫽关p,u,v,w,␾g,k,␧兴^T (9)

The numerical characteristics of the Eq.共7兲 are studied by obtaining the eigenvalues of the matrix,关⌫^⫺1(⳵E/⳵Qv)兴. The flux Jacobian⳵E/⳵Qvis given as:

(3)

A⫽

¨

u

c_␾² ␳m 0 0 共␳g⫺␳l兲u 0 0

u²

c_␾²⫹1 ␳m2u 0 0 共␳g⫺␳l兲u² 0 0 vu

c_␾² ␳mv ␳mu 0 共␳g⫺␳l兲vu 0 0 wu

c_␾² ␳mw 0 ␳mu 共␳g⫺␳l兲wu 0 0

␾gu

c_␾² ␳g␾g 0 0 ␳gu 0 0

ku

c_␾² ␳mk 0 0 共␳g⫺␳L兲ku ␳mu 0

␧u

c_␾² ␳m␧ 0 0 共␳g⫺␳L兲␧u 0 ␳mu

©

⁽¹⁰⁾

The eigenvalues of the system are derived to be:

⌳⫽共u⫹cm,u⫺cm,u,u,u,u,u兲 (11) where c_mturns out to be the well-known, harmonic expression for the speed of sound in a two-phase mixture and is given as:

1

c_m²⫽␳m

b

^␳^␾^g^c^g^g²^⫹^␳^␾^L^c^L^L²

c

⁽¹²⁾

The behavior of the two-phase speed of sound is plotted in Fig. 1 as a function of the gas porosity; at either limit the pure single- phase acoustic speed is recovered. However, away from the single-phase limits, the acoustic speed rapidly drops below either limit value and remains at the low-level in most of the mixture regime. As a consequence, the local Mach number in the interface region can be large even in low speed flows.

We re-emphasize a critical observation at this point: the equation system 共7兲–共9兲 is completely defined and does not require ad-hoc closure models for the variation of mixture density with pressure. In that respect alone this represents a significant ad- vancement over most other cavitation models in the literature. The acoustic speeds for individual phases are well-defined physical quantities, which may be specified, and so is the case with physical material densities (␳g,␳L) for each individual phase. For low pressure incompressible regimes, the material densities may be assumed to be constant without significant error in the solutions.

However, in its most general form the material densities for each phase may be obtained from the pressure using their respective physical equations of state共e.g., ideal gas law for gases, etc.兲 if that is so desired. If temperature variations were significant, this would involve solving an additional equation for the mixture en- ergy共this formulation is to be presented in a future paper兲.

To obtain an efficient time-marching numerical scheme, preconditioning is now applied to the system in Eq.共7兲, in order to rescale the eigenvalues of the system so that the acoustic speeds are of the same order of magnitude as the local convective velocities. This is achieved by replacing⌫ in Eq. 共7兲 by ⌫p.

⌫p

⳵Qv

⳵t ⫹

⳵E

⳵x ⫹

⳵F

⳵y ⫹

⳵G

⳵z ⫽S⫹Dv (13)

where

⌫p⫽

l

␤

c_␾² 0 0 0 共␳g⫺␳L兲 0 0

␤u

c_␾² ␳m 0 0 共␳g⫺␳L兲u 0 0

␤v

c_␾² 0 ␳m 0 共␳g⫺␳L兲v 0 0

␤w

c_␾² 0 0 ␳m 共␳g⫺␳L兲w 0 0

␤␾g

c_g² 0 0 0 ␳g 0 0

␤k

c_␾² 0 0 0 共␳g⫺␳L兲k ␳m 0

␤␧

c_␾² 0 0 0 共␳g⫺␳L兲␧ 0 ␳m

m

Here the parameter ␤ has been introduced to precondition the eigenvalues. The modified eigenvalues of the preconditioned system are given as:

Fig. 1 Speed-of-sound in a two-phase gas-liquid mixture

(4)

⌳P⫽

冉

^u²

冉

¹^⫹^␤¹

冊

^⫹c^m^⬘^,^u²

冉

¹^⫹¹^␤

冊

^⫺c^m^⬘^,u,u,u,u,u

冊

⁽¹⁴⁾

where

c_m⬘⫽1

2

冑

^u²

冉

¹^⫺^␤¹

冊

²^⫹4^c^␤^m² ⁽¹⁵⁾

Equation 共15兲 indicates that by setting ␤⫽(cm

2/u_p²) where u_p

⫽max(u,0.01cm) the pseudo-acoustic speed is of the order of u at all mixture composition values. We note that, a rigorous definition of the physical two-phase acoustic speed c_m 共as has been done here兲 is critical to obtaining noise-free propagation of pressure waves across interfaces where the density and acoustic speed are varying rapidly.

Cavitation Source Terms

In the present effort, the cavitation source term is simplified via a simplified nonequilibrium, finite rate form as follows:

Sg⫽Kf␳L␾L⫹Kb␳g␾g (16) where the constant K_f is the rate constant for vapor being generated from liquid in a region where the local pressure is less than the vapor pressure. Conversely, K_bis the rate constant for reconversion of vapor back to liquid in regions where the pressure exceeds the vapor pressure. Here, the rate constants are specified using the form共given by Merkle et al. 关5兴兲 as follows:

Kb⫽

b

^␶¹^b

^冉

⁰^Q^L^⬁^⬁

^冊冋

¹²^p^␳^⫺p^⬁^Q^v^⬁²

册

^p^p^⬍p^⬎p^v^v

c

(17)

Kf⫽

冋

^␶¹^f

^冉

⁰^Q^L^⬁^⬁

^冊冋

¹²^p^␳^⫺p^⬁^Q^v^⬁²

册

^p^p^⬎p^⬍p^v^v

册

p_v⫽p⬁⫺1

2␳⬁Q_⬁²* Cav. No.

␶f⫽time constant for vapor formation

␶b⫽time constant for liquid reconversion Cav. No.⫽p_⬁⫺pv

1 2␳⬁Q_⬁²

We note that the nonequilibrium time scales have not been correlated with experimental data. For steady attached cavitation this simplified form may be adequate since the cavitation time scales do not interact with the fluid time scales if the cavitation rate constants are fast enough. This point has been demonstrated by repeating a calculation for a given cavitation number at three different rates共see Results section兲. For unsteady cavitation, however, the details of how the nonequilibrium source term is specified could be crucial since it may couple with transient pressure waves. The development of a more rigorous nonequilibrium source term model is a topic of ongoing research.

Turbulence Models

The standard high Reynolds number form of the k-␧ equations form the basis for all turbulence modeling in CRUNCH. Transport equations for the turbulent kinetic energy and its dissipation rate are solved along with the basic momentum and energy equations.

These equations, with supplemental low Re terms, are as follows,

⳵␳k

⳵t ⫹

⳵

⳵xi

冉

^␳uⁱ^k^⫺

冉

^␮⫹^␮^␴^T^k

冊

^⳵x^⳵kⁱ

冊

^⫽P^k^⫺^␳␧⫹S^k

⳵␳␧

⳵t ⫹

⳵

⳵xi

冉

^␳uⁱ^␧⫺

冉

^␮⫹^␮^␴␧^T

冊

^⳵x^⳵␧ⁱ

冊

^⫽c¹^f¹^P^k^⫺C²^f²^␳␧⫹S^␧

(18) Pk⫽␶i j

⳵ui

⳵xj

, ␶i j⫽⫺2

3␳k⫹2␮T*

冉

^S^{i j}^⫺¹³^⳵u^⳵x^k^k^␦^{i j}

冊

^,

␮T⫽C␮f_␮␳k²

␧

where,␴k, ␴␧, C1, and C2are the modeling constants, and f1, f2, f_␮are empirical modeling functions to account for low Rey- nolds number共near wall兲. 共They equal unity in the high Reynolds number form.兲

Low Reynolds number effects in the near wall are accounted for by using an extension of the near wall model of So et al.关10兴.

This model has been shown to reproduce the near wall asymptotic relations for the Reynolds stress and kinetic energy accurately. In this model, the damping functions, f₁, f₂, f_␮, are defined as follows,

f₁⫽1.0⫺exp

b

^⫺

冉

^Re⁴⁰^t

冊

²

_c

f₂⫽1⫺2

9exp

冋

^⫺

冉

^Re⁶^t

冊

²

册

S_␧⫽1

4c₃␮L

冋冉

^⳵k^⳵x^1/2

冊

²^⫹

冉

^⳵k^⳵y^1/2

冊

²^⫹

冉

^⳵k^⳵z^1/2

冊

²

册

⁽¹⁹⁾

where

Ret⫽␳k²

␮L␧

f_␮⫽共1⫹4 Ret⫺3/4兲tanh

冉

^Re¹²⁵^k

冊

where

Rek⫽␳

冑

ky

␮

The constants for this model are given as follows.

c_␮⫽0.09, ␴k⫽1.4, ␴␧⫽1.4, (20) c₁⫽1.44, c2⫽1.92, c3⫽2.9556

This model has been tested for a number of standard turbulence test cases to verify that it gives correct boundary layer growth and reattachment length for recirculating flows. For cavitating flows, where large density gradients are present due to multi-fluid composition, the rigorous validation of the turbulence model will require further study; however, this is beyond the scope of the present paper.

Unstructured Crunch Code Overview

The multi-phase formulation derived in the previous section has been implemented within a three-dimensional unstructured code CRUNCH; a brief overview of the numerics is given here and we refer the reader to Hosangadi et al.,关9,11兴 and Barth 关12,13兴 for additional details. The CRUNCH code has a hybrid, multi- element unstructured framework which allows for a combination of tetrahedral, prismatic, and hexahedral cells. The grid connec- tivity is stored as an edge-based, cell-vertex data structure where a dual volume is obtained for each vertex by defining surfaces, which cut across edges coming to a node. An edge-based framework is attractive in dealing with multi-elements since dual sur-

(5)

face areas for each edge can include contributions from different element types making the inviscid flux calculation ‘‘grid transparent.’’

The integral form of the conservation equations are written for a dual control-volume around each node as follows:

⌫p

⌬QvV

⌬t ⫹

冕

⳵⍀i

F共Qv,n兲ds⫽

冕

⍀i

SdV⫹

冕

⳵⍀i

D共Qv,n兲ds (21) The inviscid flux procedure involves looping over the edge list and computing the flux at the dual face area bisecting each edge.

A Riemann problem is solved for using higher order reconstructed values of primitive variables at the dual face. Presently, a second- order linear reconstruction procedure 共following Barth 关13兴兲 is employed to obtain a higher-order scheme. For flows with strong gradients, the reconstructed higher variables need to be limited to obtain a stable TVD scheme. We note that the inviscid flux as outlined above is grid-transparent since no details of the element type are required in the flux calculation.

The Riemann problem at the dual face is computed with a Roe- averaged, flux-differenced solver which defines the flux, F_m, as:

F_m⫽1

2

b

^⌬F^⫺^F^共Q^共Q^Roe^m^m^⫺^,n^ជ^,n^ជôiôi^兲⫹F共Q^兲⫹⌬F^m^⫹^⫹^共Q^,n^ជôi^Roe^m^兲⫹^,n^ជôi^兲

c

⁽²²⁾

where nជoi is the vector area of the dual face crossing edge eជoi, and the subscript ‘‘Roe’’ denotes Roe-averaged variables. The flux differences⌬F^⫺and⌬F^⫹are given as

⌬F^⫹共QRoe

m ,nជoi兲⫽¹2共RRoe

m 共⌳⫹兩⌳兩兲RoeL_Roe^m 兲共Qm⫹⫺Qm⫺兲 (23)

⌬F^⫺共QRoe

m ,nជoi兲⫽¹2共RRoe

m 共⌳⫺兩⌳兩兲RoeL_Roe^m 兲共Qm⫹⫺Qm⫺兲 (24) Here, R and L are, respectively, the right and left eigenvector matrices, while L is the diagonal eigenvalue matrix. The Roe av- eraged velocities and scalars are defined as

yRoe⫽y^⫹

冑

␳^⫹⫹y^⫺

冑

␳^⫺

冑

␳^⫹⫹

冑

␳^⫺ , 兵^y苸共u,v,w,k,␧兲其 ⁽²⁵⁾ while, the volume fraction, ␾g, is defined by doing a simple arithmetic average:

␾g,Roe⫽¹2共␾g⫹⫹␾g⫺兲 (26) This average value of the␾gis used to define the speed of sound of the mixture at the dual face.

For efficient computation of large 3D problems a parallel framework for distributed memory systems has been implemented, along with a time-marching implicit solution procedure.

The sparse implicit matrix is derived by doing a Euler explicit linearization of the first-order flux, and a variety of iterative sparse matrix solvers, e.g., GMRES, Gauss-Seidel procedure, are avail- able in the code共see Hosangadi et al. 关9兴 for details兲. The parallel framework is implemented by partitioning the grid into sub- domains with each subdomain residing on an independent proces- sor. The message passing between processors has been implemented using MPI to provide portability across platforms.

Results and Discussion

The multi-phase system described in the previous sections has been applied extensively to both cavitating and non-cavitating problems. Since the single-phase incompressible formulation is a subset of the more general multi-phase system, the code was first validated for standard incompressible test cases in the literature 共see Ahuja et al. 关7兴, for results兲. In the present section, we will focus on the following sheet-cavitation problems which have been studied extensively in the literature: 1兲 cylindrical headform; and 2兲 the NACA 66 hydrofoil.

All the simulations presented here have been computed as steady-state calculations, using the two-equation (k-e) turbulence model with near-wall damping described earlier. The liquid to gas density ratio was specified to be 100. The baseline forward and backward cavitation rate terms␶f,␶b, are specified to be 0.001 s unless otherwise specified. Furthermore, the calculations were performed on parallel distributed memory platforms using 8 processors. We note that although the flowfields simulated here are essentially two-dimensional in character, the CRUNCH code is three-dimensional; the simulations performed here were three- dimensional with minimal resolution in the z coordinate. Hence, we do not anticipate any serious issues in computing a truly three- dimensional flowfield apart from the increased numerical cost.

Cylindrical Headform Simulations. Simulations of the Rouse and McNown 关14兴 experiments for water flowing over a hemisphere/cylinder geometry are presented. The Reynolds number per inch for this configuration is 1.36⫻10⁵. A hexahedral grid with dimensions of 221⫻113 was used for all the cavitation numbers reported. The grid was clustered both radially near the surface as well as axially at the bend where sphere meets with the cylinder. Simulations for three cavitation numbers of 0.4, 0.3, and 0.2 are presented. The grid is best suited for the highest cavitation number of 0.4 since the cavity zone is completely contained within the zone of high clustering. Hence the results for the 0.4 case will be analyzed in-depth to study the flow features in the cavity closure region, as well as other sensitivity studies such as the effect of the rate constant.

The cavitation zone for a cavitation number of 0.4 is shown in Fig. 2 for the baseline source term rates given earlier. We note that the gas/liquid interface is sharp and captured within a couple of cells is most of the flow except at the tip on the top corner of the back end where a marginal ‘‘pulling’’ or extension of the cavitation zone is observed. This is attributed to the strong shear at this location as the flow turns around the cavitation zone to enclose it. The details of the recirculation zone in the cavitation closure region are illustrated in Fig. 3 by plotting the streamlines of the flow. A sharp recirculation zone originating from the tip of the cavitation zone is observed where a re-entrant jet is observed to transport fluid back into the cavity. As we shall see from the surface pressure distribution, the re-entrant jet generates a local pressure peak at the rear of the cavity thereby keeping this cavity shape stable.

Fig. 2 Resolution of the cavitation zone interface on the nu- merical grid

(6)

The turbulent viscosity levels共nondimensionalized by the laminar viscosity兲 generated in the flowfield are shown in Figs. 4共a兲 and 4共b兲. In Fig. 4共a兲 we plot the overall global characteristics of the turbulent wake while Fig. 4共b兲 gives a close-up view of the cavitation zone itself. From Fig. 4共a兲 we observe the incoming boundary layer lifting off as the flow turns and jumps over the cavitation zone. The cavitation zone itself is observed to be very

‘‘quiet’’ with turbulence present mainly in the interface region and the shear layer above it. The plot depicts increasing turbulence intensity levels near the cavity closure/re-entrant jet region, and a fully turbulent wake ensues. We note that since the calculation was performed with a steady RANS turbulence model, the wake simulated corresponds to the time-averaged solution which

reproduces the mean mixing rate. If the calculation were performed with an unsteady LES model, an unsteady wake would result with vortices being shed off the cavity tip.

The characteristics of the flowfield at cavitation numbers 0f 0.4, 0.3, and 0.2 are compared in Fig. 5 by plotting both the cavitation zone as well as the surface pressure distribution. Overall the surface pressure profiles compare very well with experimental data.

In particular, the cavity closure region appears to be captured well; the location, the pressure gradient, and the over- pressurization magnitude in the wake are very close to the data for the two higher cavitation numbers of 0.4 and 0.3. For the 0.2 case, the cavity zone is very large and extends into the region where the grid is not as finely clustered, and we attribute the slight under Fig. 3 Flow streamlines depicting recirculation zoneÕre-entrant jet in the cavity closure region

Fig. 4 Turbulent viscosities„␮tÕ␮L…contours in the headform cavity flowfield

(7)

prediction of the cavity length to grid quality. The coarser grid in the interface region is also reflected in the broader interface zone.

It should be noted that all three cases were computed with the baseline source term rates.

To evaluate the impact of the source term rates␶f and␶b, we computed the 0.4 cavitation case using the following three different rates: 1兲 the baseline rate 共␶f⫽0.001 s, ␶b⫽0.001 s兲; 2兲 a faster rate (baseline⫻10); and, 3兲 a slower rate (baseline⫻0.1).

The surface pressure comparisons are shown in Fig. 6 for the three different rates. The results for the baseline rate and the faster rate are very close; however, as expected, the faster rate gives slightly better results in the wake with the cavity closure being sharper.

The slower rate (baseline⫻0.1) gives a more diffused cavity particularly in the rear with the time scale for the nonequilibrium source term coupling with the convective time scale. As the cavity

closure region gets more diffuse, the strength of the re-entrant jet weakens, and consequently, the pressure does not exhibit the pattern of over-pressurization and relaxation which is observed in the data as well as the other two faster rate cases. We note that the baseline rates have been used to compute all the results reported here including the hydrofoil case. At least for steady-state problems, the nonequilibrium source terms perform adequately and the precise value of the rate appears to be problem/grid independent as long as the time scale is sufficiently fast enough共as defined by our baseline rate兲 in comparison to the characteristic convective time scales.

NACA 66 Hydrofoil Simulations. In this section, we present results for sheet cavitation on a NACA 66 hydrofoil 共Shen and Dimotakis 关15兴兲. The freestream flow is at a 4 degree angle-of- Fig. 5 Cavitation zone and surface pressure profiles at various cavitation numbers for hemisphereÕcylinder headform

(8)

attack, with a Reynolds number of 2⫻10⁶based on chord length.

The hybrid grid used for the calculations is shown in Fig. 7; it has a combination of approximately 200,000 tetrahedral and prismatic cells with clustering around the surface of the hydrofoil. Calcula- tions for two different cavitation numbers of 0.84 and 0.91 are discussed in the following paragraphs. Note that the calculations shown here were computed with a wall function procedure.

The pressure contours for the flowfield at a cavitation number of 0.84 are plotted in Fig. 8. We observe that the pressure contours cluster around the cavitation boundary where the density gradient is very large and the flow turns around the cavitation bubble. The gas void fraction contours at the two cavitation num-

bers are shown in Fig. 9. At the lower cavitation number of 0.84, the cavity extends up to 50 percent of chord. As the cavitation number increases the gas bubble region decreases in length and comes closer to the surface; the void fraction contours at a cavitation number of 0.91共Fig. 9共b兲兲, indicate that the bubble extends to only 30 percent of chord. The thin confinement region of the gas bubble highlights the potential for using local grid adaption within our unstructured framework.

The surface pressure profiles at the two cavitation numbers are plotted in Figs. 10共a兲 and 10共b兲 for comparison the experimental Fig. 6 Sensitivity of cavitation solution to the cavitation

source term rate

Fig. 7 The prismaticÕtetrahedral grid used to capture cavita- tion on the NACA 66 modified hydrofoil

Fig. 8 A representative pressure distribution on the NACA 66 hydrofoil at 4 degrees angle of attack and cavitation number of 0.84

Fig. 9 Cavitation bubbles indicated by void fraction contours on the NACA 66 modified hydrofoil at cavitation numbers of„a… 0.91 and„b…0.84

(9)

data points as well as results from the numerical study by Singal et al.关16兴 are also shown. We note that the simulation by Singal et al.关16兴 was computed on a structured hexahedral grid using a user specified equation of state, and had a pdf based model for cavitation. We make the following observations; 1兲 Both the numerical simulations are very similar to each other despite the differences in the numerics and the formulation. The good comparison between the two numerical results is an important validation of the more fundamental formulation presented here which did not require a user specified closure model for the equation of state; 2兲 While the comparisons with the experimental data are reasonably good, the numerical results underpredict the length of the cavity and overpredict the aft recovery pressure. A preliminary investi- gation indicates that this may be related to a combination of local grid resolution in the interface region, as well as the details of the near wall turbulence model. However, for the purposes of illus- trating the applicability of the new formulation the current results were deemed to be acceptable.

Conclusions

A multi-phase model for low speed gas-liquid mixtures has been developed by reducing the compressible system of equations to an acoustically accurate form that performs efficiently in the incompressible regime. In particular, the equation system does not require ad-hoc density-pressure relations to close the system and yields the correct acoustic speed as a function of local mixture composition. For efficient steady-state solutions, the physical mixture acoustic speed is preconditioned to obtain good convergence.

The solution procedure has been implemented within a hybrid, multi-element unstructured framework which operates in a parallel, domain-decomposed environment for distributed memory systems.

Detailed results are presented for steady-state sheet cavitation in an hemisphere/cylinder geometry as well as a NACA 66 hydrofoil at various cavitation numbers. Good comparison is obtained with experimental data, and in particular the details of the cavity closure, and the re-entrant jet are captured well in the simulation. Examination of the turbulence characteristics indicates that turbulent kinetic energy associated with the upstream boundary layer jumps over the cavity, and the interior of the cavity itself is nearly laminar. However, the recirculation zone in the re-entrant jet region generates a fully turbulent wake which stabilizes the cavity zone and generates a pressure pattern of over-pressurization and gradual relaxation in the wake. Sensitivity of the cavitation zone to the finite rate source terms was examined. The solution was found to be relatively insensitive to the source term rate as long as the rate was high enough. As the source term rate dropped by an order of magnitude, coupling between the convective time scale and the source term was observed leading to a more diffuse cavity in the closure region and poorer pressure predictions.

Clearly, for unsteady detached cavitation problems, the modeling of the source term needs to be examined carefully and this is an area of ongoing work.

References

关1兴 Kubota, A., Kato, H., and Yamaguchi, H., 1992, ‘‘Cavity Flow Predictions Based on the Euler Equations,’’ J. Fluid Mech., 240, pp. 59–96.

关2兴 Chen, Y., and Heister, S. D., 1996, ‘‘Modeling Hydrodynamic Nonequilibrium in Cavitating Flows,’’ ASME J. Fluids Eng., 118, pp. 172–178.

关3兴 Delaunay, Y., and Kueny, J. L., 1990, ‘‘Cavity Flow Predictions based on the Euler Equations,’’ ASME Cavitation and Multi-Phase Flow Forum, Vol. 109, pp. 153–158.

关4兴 Janssens, M. E., Hulshoff, S. J., and Hoeijmakers, H. W. M., 1997, ‘‘Calcu- lation of Unsteady Attached Cavitation,’’ AIAA-97-1936, 13th AIAA CFD Conferences, Snowmass, CO, June 29–July 2.

关5兴 Merkle, C. L., Feng, J. Z. and Buelow, P. E. O., 1998, ‘‘Computational Mod- eling of the Dynamics of Sheet Cavitation,’’ Proceedings of the 3rd Interna- tional Symposium on Cavitation, Grenoble.

关6兴 Kunz, R. F., et al., 1999, ‘‘A Preconditioned Navier-Stokes Method for Two- Phase Flows with Application to Cavitation Prediction,’’ 14th AIAA CFD Conference, Norfolk, VA, June 28–July 1.

关7兴 Ahuja, V., Hosangadi, A., Ungewitter, R., and Dash, S. M., 1999, ‘‘A Hybrid Unstructured Mesh Solver for Multi-Fluid Mixtures,’’ AIAA-99-3330, 14th AIAA CFD Conference, Norfolk, VA, June 28–July 1.

关8兴 Hosangadi, A., Sinha, N., and Dash, S. M., 1997, ‘‘A Unified Hyperbolic Interface Capturing Scheme for Gas/Liquid Flows,’’ AIAA-97-2081, 13th AIAA CFD Conferences, Snowmass, CO, June 29–July 2.

关9兴 Hosangadi, A., Lee, R. A., Cavallo, P. A., Sinha, N., and York, B. J., 1998,

‘‘Hybrid, Viscous, Unstructured Mesh Solver for Propulsive Applications,’’

AIAA-98-3153, AIAA 34th JPC, Cleveland, OH, July 13–15.

关10兴 So, R. M. C., Sarkar, A., Gerodimos, G., and Zhang, J., 1997, ‘‘A Dissipation Rate Equation for Low Reynolds Number and Near-Wall Technique,’’ Theor.

Comput. Fluid Dyn., 9, pp. 47–63.

关11兴 Hosangadi, A., Lee, R. A., York, B. J., Sinha, N., and Dash, S. M., 1996,

‘‘Upwind Unstructured Scheme for Three-Dimensional Combusting Flows,’’

J. Propul. Power, 12, No. 3, May–June, pp. 494–503.

关12兴 Barth, T. J., 1991, ‘‘A 3-D Upwind Euler Solver for Unstructured Meshes,’’

AIAA-91-1548, 10th AIAA CFD Conference, Honolulu, HI, June.

关13兴 Barth, T. J., and Linton, S. W., 1995, ‘‘An Unstructured Mesh Newton Solu- Fig. 10 Surface pressure distribution on the NACA 66 hydro-

foil using wall-function procedure

(10)

tion for Compressible Fluid Flow and Its Parallel Implementation,’’ AIAA-95- 0221, 33th AIAA Aerospace Sciences Meeting at Reno, NV, Jan 9–14.

关14兴 Rouse, H., and McNown, J. S., 1948, ‘‘Cavitation and Pressure Distribution:

Head Forms at a Zero Angle of Yaw,’’ Technical Report: State University of Iowa Engineering Bulletin No. 32.

关15兴 Shen, Y. T., and Dimotakis, P., 1989, ‘‘The Influence of Surface Cavitation on

Hydrodynamic Forces,’’ 22nd American Towing Tank Conference, St Johns, NF, August 8–11.

关16兴 Singal, A. K., Vaidya, N., and Leonard, A. D., 1997, ‘‘Multi-Dimensional Simulation of Cavitating Flows Using a PDF Model for Phase Change,’’

FEDSM97-3272, 1997 ASME Fluids Engineering Division Summer Meeting, Vancouver, British Columbia, Canada, June 22–26.

Simulations of Cavitating Flows Using Hybrid Unstructured

Vineet Ahuja

Ashvin Hosangadi

Srinivasan Arunajatesan

Simulations of Cavitating Flows Using Hybrid Unstructured

Meshes

冉

冊

冉

冊

冉

冊

冉

冊

¨

©

¨

©

b

c

l

m

冉

冉

冊

冉

冊

冊

冑

冉

冊

b

冉

冊 冋

册

c

冋

冉

冊 冋

册

册

冉

冉

冊

冊

冉

冉

冊

冊

冉

冊

b

冉

冊

c

冋

冉

冊

册

冋 冉

冊

冉

冊

冉

冊

册

冉

冊

冑

冕

冕

冕

b

c

冑

冑

冑

冑

^冉

^冊冋

^冉

^冊冋

_c

冋冉