Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
24th Symposium on Naval Hydrodynamics Fukuoka, JAPAN, 8-13 July 2002 Prediction of Vortex Cavitation Inception Using Coupled Spherical and Non-Spherical Models and UnRANS Computations C.-T. Hsiao and G. L. Chahine (DYNAFLOW, INC., USA) ABSTRACT A spherical and a non-spherical bubble dynamics model were developed to study cavitation inception, scaling, and dynamics in a vortex flow. The spherical model is a modified Rayleigh-Plesset model modified to account for bubble slip velocity and for non-uniform pressures around the bubble. The non-spherical model is embedded in an Unsteady Reynolds-Averaged Navier-Stokes code with appropriate free surface boundary conditions and a moving Chimera grid scheme around the bubble. The effect of non-spherical deformation and bubble/flow interaction on bubble dynamics is illustrated by comparing spherical and non-spherical models. It is shown that non-spherical deformations and bubble/flow interactions are important for accurate prediction of cavitation inception. The surface- averaged pressure modified Rayleigh-Plesset scheme is a significantly improvement over the conventional spherical model and is able to capture the volume changes of the bubble during its capture. It presents a fast scheme to study scaling. Scaling effects on cavitation inception is preliminarily studied using two different Reynolds number due to two different chord lengths. The nuclei size effect on the prediction of cavitation inception is also studied and its important effect highlighted. I. Intro(luction Prediction and scaling of vortex cavitation on propellers from model tests is of critical importance to naval applications since this type of cavitation often occurs first and affects the discretion of the propeller. In order to do so it is essential to be able to understand and model accurately the mechanisms that enter into play in the cavitation inception process of a vortex flow. An analysis of the problem indicates that, in order to accurately predict cavitation inception, one needs to consider the bubble dynamics, the detailed viscous flow structure in the low pressure regions of the propeller flow field, and their interaction. Such an approach may lead to the conclusion that some of the effects could be neglected, but is needed before approximations can be made satisfactorily. Numerical predictions (Chahine 1990, 1994, 1995) and recent experimental observations (Arndt and Maines 2000), indicates that the cavitation inception starts with a spherical nucleus drawn into the vortex core from the free stream and which grows significantly. As the bubble approaches the vortex core, non-spherical deformations become significant. The bubble deforms, elongates, and later breaks up into multiple bubbles as it travels downstream and reaches the vortex axis. Numerical simulations of bubble capture by a tip vortex (Ligneul and Latorre 1989, 1993, Hsiao and Pauley 1999, Hsiao et al. 2000, 2002) were primarily limited to a spherical bubble model which assumes that the bubble does not interact with the underlying flow field and sees a uniform pressure around it given by the value of the liquid pressure at the location of its center (in its absence) during its capture. Since the traveling bubble is usually very small relative to the vortex core size before cavitation, the modification of the vortex flow may be neglected. The assumption of a spherical model, however, is not always appropriate because the bubble may grow to a size large enough such that non-spherical deformations occur when the bubble approaches the vortex axis. A non-spherical bubble dynamics model based on the boundary element method was first developed by Chahine (1990, 1994, 1995) to study bubble deformations during capture but neglected the effect the bubble has on the underlying vortex flow field. To further account for the bubble/vortex interaction, Hsiao and Chahine (2001) developed a non-spherical bubble dynamics model with full viscous interaction between the bubble and the vortex based on unsteady Navier-
Stokes computations and on a moving Chimera grid scheme. The bubble dynamics model has been successfully validated for known cases and applied to study bubble dynamics in a Rankine line vortex. Their results showed that non-spherical bubble deformations and bubble and flow interaction can significantly influence bubble dynamics in the vortex flow below cavitation inception. However, the unsteady Navier-Stokes computations are time- consuming and not practical for simulating the whole bubble dynamics process and accounting for the nucleus motion from upstream until drawn into the vortex core. Recognizing this fact we apply the spherical model to simulate the bubble capture and then turn on the non-spherical model when the bubble size exceeds a preset limit value to best take advantage of each model. In the current study, our model is applied to predict tip vortex cavitation inception for a canonical problem in which the tip vortex flow is generated by a finite-span elliptic hydrofoil. The flow field is first obtained by a steady-state Navier-Stokes computation and provides the velocity and pressure fields for the spherical model. To study the tip vortex cavitation inception we assume that a bubble nucleus is released upstream. The spherical model is first used to track the bubble during its capture by the tip vortex. When the non-spherical model is turn on, the flow field due to the spherical bubble motion and volume change is superimposed to the global steady-state blade flow field solution to provide initial conditions for the unsteady computations. 2. Numerical Methods 2.1 Navier-Stokes Computations To best describe the tip vortex flow field around a finite-span hydrofoil, the Reynolds- Averaged Navier-Stokes (RANS) equations with a turbulence model are solved because the RANS computations have been shown to be successful in addressing tip vortex flow (Hsiao and Pauley 1998, 1999) and general propeller, propulsor, ship, and free surface flows (Taylor et al. 1998, Wilson, et al. 19981. The three-dimensional unsteady Reynolds- Averaged incompressible continuity and Navier- Stokes equations written in non-dimensional form and Cartesian tensor notations are given as ~=0, (1) axi an +ujaxj=-axi+Re axj (2) where ui = (u, v, w) are the Cartesian components of the velocity, xi =(x,y,z) are the Cartesian coordinates, p is the pressure, Re = pu*L*/,u is the Reynolds number, u* and L* are the characteristic velocity and length, p is the liquid density, and ,u is its dynamic viscosity. The effective stress tensor Fit . . is given by: ij =~¢axj + axi a- 3 it axkk -uiuj , (3) where by is the Kronecker delta and u~u' is the Reynolds stress tensor resulting from the Reynolds averaging scheme. To numerically simulate the tip vortex flow filed around a finite-span hydrofoil, a body-fitted curvilinear grid is generated. Equations (1) and (2) are transformed into a general curvilinear coordinate system. The transformation provides a computational domain that is better suited for applying the spatial differencing scheme and the boundary conditions. To solve the transformed equations, we use and modify the three-dimensional incompressible Navier-Stokes flow solver, DF_UNCLE, derived from the code UNCLE developed at Mississippi State University. UNCLE was modified at DYNAFLOW to include bubble, cavities and free surface effects (Hsiao and Chahine, 2001, Chahine and Hsiao, 20021. The DF_UNCLE code is based on the artificial compressibility method (Chorin 1967) in which a time derivative of the pressure multiplied by an artificial-compressibility factor is added to the continuity equation. As a consequence, a hyperbolic system of equations is formed and can be solved using a time marching scheme. This method can be marched in pseudo-time to reach a steady-state solution. To obtain a time-dependent solution, a Newton iterative procedure is performed at each physical time step in order to satisfy the continuity equation. The numerical scheme in DF_UNCLE uses a finite volume formulation. First-order Euler implicit differencing is applied to the time derivatives. The spatial differencing of the convective terms uses the flux-difference splitting scheme based on Roe's method (Roe, 1981) and van Leer's MUSCL method (van Leer, 1979) for obtaining the first-order and the third-order fluxes respectively. A second-order central differencing is used for the viscous terms which are simplified using the thin-layer approximation. The flux Jacobians required in an implicit scheme are obtained numerically. The resulting system of algebraic equations is solved
using the Discretized Newton Relaxation method (Vanden and Whitfield, 1993) in which symmetric block Gauss-Seidel sub-iterations are performed before the solution is updated at each Newton interaction. DF_UNCLE is also accompanied by a Baldwin-Lomax algebraic turbulence model to model the Reynolds stresses in Equation (31. 2.2 Spherical Model When one assumes that the bubble remains spherical during its volume variations, the bubble dynamics can be described by the Rayleigh-Plesset equation (Plesset 19481. The classical Rayleigh- Plesset equation can be derived from the continuity and momentum equations of a ID radial viscous flow. In the current study, this equation is modified to account for a slip velocity between the bubble and the host liquid, and to account for the non-uniform pressure field along the bubble surface. The bubble is passively convected through the known "basic" flow field of the blade or propeller. The bubble's travel velocity is generally different from the liquid velocity due to the bubble dynamics and to the presence of local flow field pressure gradients. This results in an added pressure term due to the velocity difference between the liquid and the bubble. This term is obtained by adding the velocity potential due to the flow past a sphere to the velocity potential due to a spherically oscillating bubble and using this in the pressure balance equation. As a result, an additional term of slip velocity, (u-ub)214, where u is the liquid convection velocity and ub is the bubble travel velocity, is added to the classical Rayleigh- Plesset equation as: RR+3R2 1 2 p PA + Pgo ( R ) PenCounter By 4,uR R R (uus ~ (4) R is the time varying bubble radius, Ro is the initial or reference bubble radius, ~ is the surface tension parameter, pgo is the initial or reference gas pressure inside the bubble, and Pv is the vapor pressure. The time varying gas pressure inside the bubble, pp. is related to Ego through a polytropic compression law: pgR = P9O~ ~ (5) with the polytropic gas constant k. k = 1 for isothermal behavior is used in the present study. Pencounter is the ambient pressure "felt" by the bubble during its travel. In the classical Rayleigh-plessset equation, Pencounter is the pressure that exists in the liquid, in absence of the bubble, at the location of the bubble center. This definition is quite acceptable if the pressure field does not vary significantly around the bubble. In the case of bubble capture in a vortex, such a definition used by previous investigators, leads to unlimited bubble growth once the bubble reaches a constant intensity line vortex axis. Hsiao et al. (2000, 2002) introduced the Surface Averaged Pressure (SAP) Rayleigh Plesset equation, in which Pencounter is defined as the average of all liquid pressures at the various parts of the bubble surface. This concept enabled more realistic modeling of the bubble dynamics in line vortices. We will use here the SAP approach and compare it with the fully coupled 3D bubble dynamics. To describe the bubble trajectory during capture by the tip vortex, the motion equation described by Johnson and Hsieh (1966) is used: dt p 4 CD (Uub ) | uUb| + 3 (u rub ) R. (6) where the drag coefficient CD is given by an empirical equation such as that of Haberman and Morton (19531: 24 (1+0 197ROb63 +2.6xlO Reb ) ~ (7) eb where the bubble Reynolds number is defined as Reb = ~ (8) To determine the bubble motion and its volume variation, a Runge-Kutta fourth-order scheme is used to integrate Equations (4) and (6) through time. The liquid velocity and pressures ~ Pencounter ~ are obtained directly from the RANS computations. The numerical solution of the RANS equations, however, offers the solution directly only at the grid points. To obtain the values for any specified location (x,y, z) on the bubble we need to interpolate from the background grid. To do so, an interpolation stencil (8)
and interpolation coefficients at any specified location need to be determined at each time step. To determine the interpolation stencil from the background grid, we implemented a three- dimensional point-locating scheme based on the fact that the coordinates (x,y,z~of the bubble location are uniquely represented relative to the eight corner points of the background grid stencil: 8 8 8 x = 2, NiXi, y = ~ Ni Yi ~ Z = 2, NiZi ~ (9) i=1 i=1 i=1 where N1 = (1 - V)~1 -~1 - By, N2 = ¢~1 -~1 -are ~ 10y N3 =~1-~1-~), N4 =~1 Hi, Ns = (1 - ¢~1 -yip, N6 = ¢~1 -yip, N7 = (1-oJyr(p, N8 = YELP ¢,yr,~ are the interpolation coefficients, and (xi,yi,zi ) are the coordinates of the eight corner points of a grid stencil in the background grid. Equation (9) is solved using a Newton-Raphson method. For a bubble point to be inside the grid stencil requires that the corresponding ¢,yt,~ satisfy 0 < ~ < 1, 0 < jr < 1, 0 < ~ < 1. Once the interpolation stencil and interpolation coefficients are determined, the pressure and velocities can be obtained by using Equation (9~. 2.3 Non-Spherical Model To fully account for the interaction between the bubble and the flow filed and for non-spherical deformations, a non-spherical bubble model embedded in the UNRANS computations with appropriate free surface boundary conditions and a moving Chimera grid scheme has been developed by Hsiao and Chahine (2001~. 2.3.1 Free Surface Boundary Condition To best describe the bubble surface behavior, a general free surface boundary condition satisfying both kinematic and dynamic boundary conditions is applied. The kinematic boundary condition is the Lagrangian condition that a particle on the surface must remain on the surface. For a free surface of equation Ff A, t) = 0, this can be written: DF Dt -=n ( The general free surface dynamic boundary condition is the condition of zero shear stress and balance of normal stresses at the bubble liquid interface. Here, the stress due to the gas inside the bubble is neglected. With the same simplifications used by Batchelor (1967) for deriving the dynamic boundary condition in the Cartesian coordinate system, Hodges et al. (1996) derived a dynamic boundary condition in a curvilinear coordinate system by requiring the grid to be normal to the boundary. Following their work the current study implements the dynamic boundary condition in non- dimensional form as OU ~ ~~ OW ~2 OW a,; =-g33~g a,; +g a at _g ( 22OW+ l23W) (13) We ~ 14) with gij=d~i obj. gij= axi axj (15) axk axk ask ask (12) P Pgv Re a,; Off ~e where e is the curvature, We =pu*2 L*ly is the Weber number, ~ is the surface tension, and Pgv = (Pg +Pv - poO)/pu*2 . (16) To determine the gas pressure we assume that the amount of non-condensable gas inside the bubble remains constant and that the gas satisfies the polytropic relation pg~k=constant~where ~ is the bubble volume. 2.3.2 Moving Chimera Grid Scheme The Chimera grid scheme is a grid embedding technique which provides a conceptually simple method for domain decomposition. In this approach, structured sub-grids generated around each component in the flow field or over sub- domains of complex geometries are put together without requiring the mesh boundaries to join with the global grid in any special way. In the present study, a body-fitted sub-grid is created around the bubble and overset with a global grid which is generated for the hydrofoil. Figure 1 illustrates the details of a Chimera grid system in a two- dimensional domain. The global grid consists of three different types of points: "regular" points, "overlap" points and "hole" points. The sub-grid consists of two different types of points; "regular" points and "overlap" points. The Navier-Stokes equations are solved separately for the global grid and for each sub-grid at all "regular" points. The communication between these two grids is made by
interpolating flow variables at the "overlap" points. Data from the global grid are used in the interpolation to supply outer boundary conditions to the sub-grid. Hence the effect of the main flow field is properly imposed on the bubble. For the global grid solution to account for the bubble, points in the global grid are blanked out to form a "hole" within some neighborhood of the bubble. On the fringes of this blanked-out region, data from the sub-grid solution are used for interpolation to provide interior boundary conditions for the global grid. This way the presence of the bubble and its effect are transferred to the main flow field. The implementation of the Chimera grid scheme includes two tasks: 1) identify the "hole" points and the "overlap" points, and 2) determine the interpolation stencil and interpolation coefficients. The "hole" is here defined by a "hole boundary" which consists of the grid surface at a constant value of (. A global grid point is considered to be a "hole" point if it is inside the "hole boundary". A grid point is considered to be inside the "hole" if the dot product between the vector from the closest point at the "hole boundary" to the grid point and the normal vector at the "hole boundary" at the closest point is negative or zero. The overlap points in the global grid are then determined based on the fringes of the "hole". Here, two layers of grid stencil are identified as the "overlap" points in order to enable implementation of a third order boundary condition. The method for determining the interpolation stencil and the interpolation coefficients from the background grid for each overlap point is the same as that described in section 1.2. In addition to the two tasks described above, a blanking technique is implemented in the Navier-Stokes solver to account for "holes" and interior boundaries and to account for boundary data transferred from interpolated data sets rather than from the usual boundary condition routines. In the iterative solution routines, a variable ibLi,j,k) is introduced to blank out the "hole" and "overlap" points such that no update of the variables takes place at these points: ib~i k'_J1, if apointis not blanked <16> 'J' )0, if a point is blanked In DF_UNCLE, a system of linear equation is written to solve for the flow variable difference, /` Q. between iterations and to update the flow variables by Qn+~=Qn+AQ. To ensure that AQ is zero at blanked points, we multiply the matrix coefficients of the left hand side and the vector of right hand by ib end assign the diagonal terms of the matrix to be one for Figure 1. Illustration in 2D of a Chimera grid system. Localization of the "hole" points and of the "overlap " points marked as · . 2.3.3 Initial conditions It is important to specify appropriate initial conditions for the unsteady Navier-Stokes computations when the non-spherical model is turn
on. In the current study, the initial condition is constructed by superimposing the local perturbation flow due to the bubble dynamics and the global tip vortex flow obtained by the steady-state Navier- Stokes computations. The combined velocity field, Hi =(u,v,w), for the initial condition is, therefore, obtained by R2 Ui=Ui+ 3 R(Xi-Xio) i=1,3, (17) where r is the distance from the bubble center x~O = (xO'yO,zO) to the field point xi. The combined pressure field, p, is prescribed as: R( ~ ~ RR2 [1 (R )3~ (18) where pO is the liquid pressure at the bubble center and _1 PR = PgV - R R - W ' (19) is the liquid pressure at the bubble wall. It is noted that the variables in Equations (17~-~19) are all non- dimensional. 2.4 Computational Domain and Grid Generation To compute the flow around the finite-span elliptic hydrofoil we generate an H-H type grid with 12 blocks for a computational domain which has all far-field boundaries located six (6) chord lengths away from the hydrofoil surface. For the unsteady computations, we generate an O type grid on the bubble surface (41 x 21~. Also, an O-O type 3D domain grid ~ 41 x 21 x 25 ~ is generated to extend the grid from the bubble surface to the outer boundary that we locate at25R,~s. Rns is the bubble radius when the non-spherical model is turn on. Since the bubble moves and its surface deforms during the numerical simulations, an efficient grid generation scheme that can automatically generate an appropriate grid based on the moving boundary at each time step must be integrated with the Navier- Stokes solver. We apply a grid generation scheme combining both algebraic and elliptic grid generation techniques as described by Hsiao and Pauley (19961. This grid generation scheme is selected because the algebraic grid generation technique is suited to create grid clustering and boundary-orthogonal grids at the bubble surface. This is important for resolving the flow field near the bubble surface and for applying appropriate free surface boundary conditions. The elliptic grid generation technique is then applied to smooth out the rough algebraic grid. With the current grid generation, the first grid spacing from the bubble surface is kept as a constant distance (= 0.05Rns) during the computations. To avoid uneven distribution of the surface grid points due to bubble deformation, a re-gridding scheme is applied every few time steps at the bubble surface to restore the grid distribution Cording to the initial grid-stretching factor. 3. Results and Discussion 3.1 Single Phase Flow To illustrate the method we consider an elliptic hydrofoil having a NACA 16020 cross- section with an aspect ratio equal to 3 (based on semi-span). For the steady-state computations, a 12- block grid with a total of 2.7 million grid points was generated. To optimize the grid resolution for the tip vortex, the grid clustering is adjusted using repeated computations to align the clustering around the location of the tip vortex centerline. The final refined grid has at least 16 grid points in the spanwise direction and 26 grid points in the crosswise direction within the vortex core. This results in the smallest grid cell in the vortex center having a size of about 2xlO~ CO, where CO is the hydrofoil chord length at the base section. The streamwise grid resolution near the tip region is kept almost constant end of a spacing of 2xlO 3 CO. The first grid in the normal direction to the foil surface is located 1 x 10 s CO away from the solid surface. NACA16020 oc=1 2° ~' . ' :// Re=1A4x10' ;' ~ Re=2.88x106 _ .' 0.1 0.15 0.2 0.25 X Figure 2. Pressure coefficient variations along the NACA16020 elliptic foilfor two values of the Reynolds number.
The flow field at an angle of attack equal to 12° was computed for various Reynolds numbers. The computed pressure coefficients along the tip vortex centerline for the two Reynolds numbers, Re =1.44xlO6and Re =2.88x106 are shown in Figure 2. It is seen that the locations of the minimum pressure for both these cases are very close to the hydrofoil tip at x/CO=0.05andO.04 respectively. The corresponding minima m pressure coefficients are -Cpmin =2.27 and 2.73. It is a common, even though not always accurate (Hsiao et al. 2000, 2002) to equate the cavitation inception number with the negative value of Cpmin. The ratio between the two values of Cpmin is usually presented as the ratio of the two Reynolds numbers to some power ~ usually taken to be 0.4. There are recent indications that this power should itself vary with the Reynolds number decaying at the larger values (Shen et al. 2001a). Here we see that the ratio corresponds to a value of a =0.27. Further improvements including grid refinement and higher order turbulent modeling is needed before we draw conclusions for that matter. 3.2 Bubble Dynamics with the Spherical Model For computations of bubble dynamics using the spherical model, all water properties are defined at 20° C and nuclei present in the liquid upstream are convected in the tip vortex flow. As described in the previous studies of Hs~ao and Pauley (1999) and Manes and Arndt (1999), only nuclei which pass through a restricted space volume ("window of opportunity") can be drown into the vortex center and grow explosively. This "window of opportunity" can be determined by releasing nuclei upstream and tracking their trajectory to see if they encounter the minimum pressure in the tip vortex. In the spherical computations shown here all bubble were released fromx/ CO = -0.1. The release locations in the y and z directions were chosen such that the bubbles encounter the minimum pressure in the tip vortex. Figure 3 shows a typical spherical bubble trajectory. The figure also illustrates the bubble size at the various locations along its trajectory. Figure 4 shows the bubble behavior during its capture by the tip vortex using the SAP spherical model. The case shown is for Rip =50pm, c; =2.5 and Re = 2.88x106 based on a flow speed, up = 2.88 m/s and a chord length, CO =1 m. The pressure encountered by the bubble, Pe72Counter' is shown as a function of time. Also shown is the history of the bubble radius and of the acoustic pressure generated by the bubble at the location x/Co=O,y/CO=Qz/CO=0.3. It is seen that the bubble grows continuously during its capture and reaches its maximum size after it passes the minimum pressure location. The bubble then starts to collapse as it sees an increasing encounter pressure and executes strong volume oscillations which lead to a strong acoustic signal Definition of the cavitation inception can be based on either the bubble size exceeding a limit value or on the emitted acoustic pressure exceeding a threshold. We will apply the former definition in the present study. ~ Cavi~ng bubble: i: ~::~: ~ :~ ... . .... ~ Figure 3. Bubble trajectory during its capture and resulting bubble size from actual computation 0.002s 0.002 ,~.0015 ~ 5\ 0.001 0.0005 Re=2.88x1 o6 RO=50pm c,=2~ 0.04 0.06 0.08 0.1 Time (see) 13000 12000 1 1 000 1 0 0002 _ 9000 ~ a) 8000 cut 7000 Is 6000 5000 4000 U] 3000 2000 Figure 4. Example computation of bubble dynamics for encountered pressure, bubble radius and emitted acoustic pressure versus time during its capture in the tip vortex of a NACA 16020 foil using the SAP spherical model.
3.3 Bubble Dynamics for Non-Spherical Model To include non-spherical bubble deformations and bubble/vortex interaction effects on bubble dynamics, the non-spherical bubble model is turn on when the bubble size exceeds a preset value. For consistency with the numerical scheme, this limit size is set to correspond to the smallest grid size in the tip vortex core region (Rns=2xlO 4Co). To reduce the unsteady Navier-Stokes computational burden, a computational sub-domain which contains 51x41x51 grid points is extracted from the overall computational domain as shown in Figure ~ This sub-domain is chosen large enough so that the flow quantities at the boundary of the domain are not influenced by the bubble dynamics. Figure 5. A view of the sub-domain used for the unsteady RANS computations with bubble deformation. An so type grid is then created for the bubble overset grid for the unsteady Navier-Stokes computations. To start these computations, the velocity and pressure field obtained from Equations (17) and (18) are applied as initial conditions. The initial values at all boundaries of the sub-domain are maintained over time during unsteady computations except for the solid boundary at which no-slip condition is enforced. Figure 6 shows the initial pressure contours on three selected grid planes (two of them for the global grid aid one for the overset grid) with the overlap points blanked out. It is seen that the flow field is only locally altered by the presence of the bubble. Figure 6. Init~alpressure contours on three selected grid planes (two of them for the global gad and one for the overset grid) with the overlap points blanked out. Figure 7 illustrates the bubble shapes several time steps obtained for Ro = 50pm, cr = 2.5 and Re =2.88x106 while the bubble is traveling along the tip vortex. It is seen that the bubble elongates in the axial direction, then as in the spherical model, the bubble starts to collapse after reaching its maximum size. The unsteady 3D computations, however, fail so far to continue once strong deformations develop over the bubble surface during the collapse. Comparisons of the bubble radius versus time for the spherical models (the conventional and the SAP model) and the non- spherical model are shown in Figure 8. For the non- spherical model we have used the equivalent bubble radius based on the bubble volume. It is seen that the conventional spherical model very significantly over-
predicts the maximum bubble size since the pressure, P. applied in that spherical model is the pressure existing at the bubble center. Such a model obviously does not account for pressure variations around the bubble surface. Our surface averaged pressure (SAP) spherical modelenables a more realistic evaluation of the bubble dynamics. Here, the average is obtained using the six polar points at the bubble surface, the solution of our spherical model agrees very well with the nonspherical model as shown in Figure 8. 1 tax my:: Figure 7. Bubble shapes at several time steps obtainedfor Ro = 50 pm, ~ = 2.5 arid Re =2.88x106. R=2.88xl o6 RO=50.um cs=2.5 n nn7 n nnF 0.005 0.004 0.003 0.002 0.001 Sphertcal Bubble NoSAP - S phedcal Bubble With SAP Non-Sphelical Model I .... I.... ~ .A . 0.0 4 0. 05 0 .06 0.07 0.0 8 0 .09 0.1 Tune (see) Figure 8. Comparisons of the bubble radius versus time for the spherical models (the conventional and the SAP model) and the 3D UnRANS computations. 3.4 Scaling Effect on Cavitation Inception To study scaling effects on cavitation inception, the bubble dynamics models were applied to predict cavitation inception in two different scales,C0=lm and 0.5m. All computations were conducted with the same free stream velocity, u<= = 2.88 m/s, and initial bubble radius, Ro = 50pm . As mentioned in Section 3.2, the bubble size is used for defining cavitation inception. Here, the bubble size is represented by the bubble radius for the spherical model and by the equivalent radius for the non-spherical model. A large series of computations was conducted with the spherical model (since it is very fast to run) for a range of cavitation number and is shown in Figures 9 and 10. For selected cavitation numbers computations were also conducted with the 3D UnRANS non-spherical model. Both conventional and SAP spherical models were used. Figure 9 and 10 shows the maximum size the bubble reaches at each cavitation number studied for Re = 1 .44x 1 o6 and Re = 2.88 x 1 o6 respectively. It is seen that in each case the bubble experiences an explosive growth when the cavitation number is smaller than a limit cavitation inception value. The predictions of the 3D non-spherical model (at least in terms of the equivalent radius) are very close to those obtained by the SAP spherical model. Averaging the encounter pressures around the bubble surface appears to be a justified approach. Close to the cavitation inception number, a sudden change in the maximum bubble size is seen, in which case direct punctual comparison between the three models could give the impression of very different results. To determine the cavitation inception from the current result, we must define the cavitation inception event first. If the cavitation inception event is defined based on the sudden change in the maximum bubble size, the cavitation inception numbers for both scales are read as 2.13 and 2.58 respectively. This indicates that the cavitation inception number varies between two scales as the ratio of the Reynolds numbers to the same power as the ratio of the two Cpmin. In other words, ~2/~1 ~ CPmin2/CPmin1 = (Re2/Re1 )0 28 is similar to the prediction of single phase flow. However, if the cavitation inception event is based on when the maximum bubble size reaches some value, say 2mm, then the prediction of the cavitation inception number may be quite different, e.g. Hi ~ Re036 if the criterion for cavitation inception event is defined for RmaX > 2mm.
3.5 Nuclei Size Effect on Cavitation Inception The nuclei size is already known to have a profound effect on the prediction of tip vortex cavitation inception (Hsiao and Pauley, 1999 and Shen et al., 2001b. Hsiao et al. 2000, 2002~. Previous numerical studies, however, relied on a spherical model. Here, nuclei size effect is also studied using our non-spher~cal model. 3D computations were conducted for two different nuclei sizes, Ro = 20pm arid 50,um, at the same Reynolds number, Re = 2.~8x106 . Comparison between the bubble behavior with these two initial nuclei sizes are also shown using the curves maximum bubble size versus cavitation number. It is seen from Figure 11 that the smaller nuclei requires a much smaller value of the cavitation number for cavitation inception. Then, for values of cavitation number lower than the inception limit, a more abrupt change in the RmaX versus s curve is seen. However, as already known, we find again that the maximum bubble size reached by the bubble is independent of the initial nuclei size for cavitation number much smaller than the limit cavitation inception value. As we did for Rip =50,um, we compare in Figure 11 the SAP spherical model computations and the 3D non-spherical computations for Rid = 20pm. It is seen that the SAP spherical model over-predicts a little the maximum bubble growth size at low cavitation inception numbers giving larger bubble sizes. This may reflect that using only 6 polar points may not be very accurate for large deformations and bubble sizes. 4. Conclusions A bubble dynamics model combining a spherical model and a non-spherical model embedded in UnRANS computations was developed to predict cavitation inception in complex flow configurations and was applied here for tip vortex cavitation inception on a finite-span hydrofoil. Comparisons between the classical spherical model and our current model showed that the non- spherical deformation and bubble/flow interaction can be very important. A Surface Averaged Pressure (SAP) scheme can be applied to improve significantly the spherical model when predicting the cavitation inception. Scaling effects were demonstrated by comparing the results for two different scales. It is again found as in Hsiao et al. (2000, 2002) that the definition of cavitation inception or the means for detecting inception, can significantly affect the predicted cavitation inception number. Comparison between two different nuclei sizes showed that smaller nuclei sizes result in smaller inception values, and that the maximum size at low 6 is independent of the initial nuclei size. 0.01 ~ nmn 0.008 ~ ~7 o06 .005 .004 0.003 ao02 0.001 Re=1 A4x106 RO=5Olim CO=O~m ~ \ Spherical M odel N 0 SAP - \ Spherical M odel Wile SAP ~ Non-Spheri:al M odel O -I ~ ~ , I . ~ ~ ~ I A___. . . 1 .9 2 2.1 2.2 23 2.4 Cavitation no. Figure 9. Comparison of the maximum bubble radius size versus cavitation number between spherical (conventional and SAP) and non- spherical models for the small scale, CO = 0.5m, Re = 1.44 xl 06, Rg = 50,um. 0.0125 0.01 E 0.0075 c: 0.005 0.0025 Re=2.88:~106 RO=50,um CO=lm \ ~ -- Spherical Mod b No SAP \ Sphedcal Mod b Wffl SAP ~~ Non Spherical Model 2.3 2.4 2.5 2.6 Cavitation no. 27 2.8 Figure 10. Comparison of the maximum bubble radius size versus cavitation number between spherical (conventional and SAP) and non- spherical models for the large scale, CO = I m, Re = 2.8 8x Jo6 ~ Ro = 50,um.
Re=2. 88x1 o6 Cn=1 m 0.005~' 0.004 X 0.003 a: 0.002 0.001 ~ - - R. =20 Lun SDhertal ModelWlthSAP _ - _ ~ _ . ~ I · · __~ 2.2 2.3 2.4 2.5 2.6 2.7 Cavitation No. u Ro-50 rum Spherical Model With S AP Ro_20 run Nm-SphericalMod~ \ ~ Ro=SO run Na~-Spherical Modd '\ \ Figure 11. Comparison of the maximum bubble radius size versus cavitation number between spherical (conventional and SAP) and non- spherical models for two different nuclei sizes, Ro =20 arm 50pm at Re = 2.88x 106. ACKNOWLEDGMNETS This work was conducted at DYNAFLOW, INC. (www.dvnaflow-inc.com) and has been supported by the Office of Naval Research under contract No. N0014-99-C-0369 monitored by Dr. Ki- Han Kim. This support is greatly appreciated. REFERENCES t1] Arndt, R.E and Maines, "Nucleation and Bubble Dynamics in Vortex Flows," ASME Journal of Fluids Engineering, Vol. 122, 2000, pp. 48 493. [2] Batchelor, G. K., Fluid Dynamics, Cambridge University Press, 1967. [3] Chahine, G.L., "Nonspherical Bubble Dynamics in a Line Vortex," ASME Cavitation and Multiphase Flow Forun~ Toronto, Canada, Vol. 98,Junel990,pp.121-126. [4] Chahine, G.L, "Bubble Dynamics and Cavitation Inception in Non-Uniform Flow Fields" 20th Svmposium on Naval Hvdrodvnamics, Santa Barbara, California, August 1996, pp. 290-311. t5] Chahine, G. L., "Bubble Interaction with Vortices," Vortex Flow, Chapter 18, Ed. S. Green, Kluwer Academic, 1995. [6] Chahine, G.L., Hsiao, C.-T., " Numerical Investigation of Sheet and Cloud Cavitation Inception and Dynamics on Propeller Blades," Rep. 98003- 1, 2002, Dynaflow, Inc., Jessup, MD. [7J Chorin, A. J., "A Numerical Method for Solving Incompressible Viscous Flow Problems," Journal of Computational Physics, Vol. 2, 1967, pp. 12-26. t8] Haberman, W.L., Morton, R.K., "An Experimental Investigation of the Drag and Shape of Air Bubbles Rising in Various Liquids," Report 802, 1953, DTMB. [9] Hodges, B., Street, R. Zang, Y., "A Method for Simulation of Viscous, Nonlinear, Free-Surface Flows," 20th Svmposium on Naval Hydrodynamics, 1996,pp. 791-809. t103 Hsiao, C.-T., "Numerical Study of the Tip Vortex Flow Over a Finite-Span Hydrofoil, " Ph.D. Thesis, Department of Mechanical Engineering, The Pennsylvania State University, Adviser L.L. Pauley, 1996 [11] Hsiao, C.-T., Pauley, L.L., "Numerical Study of the Steady-State Tip Vortex Flow over a Finite- Span Hydrofoil," ASME Journal of Fluid Engineering, Vol. 120, 1998, pp. 345-349. [12] Hsiao, C.-T., Pauley, L.L., "Numerical calculation of Tip Vortex Flow Generated by a Marine Propeller," ASME Journal of Fluid Engineering, Vol. 121, 1999a, pp. 63~645. t13] Hsiao, C.-T., Pauley, L.L., "Study of Tip Vortex Cavitation Inception Using Navier-Stokes Computation and Bubble Dynamics model," ASME Journal of Fluid Engineering, Vol. 121, l999b, pp. 198-204. [14] Hsiao, C.-T., Chahine, G.L., Liu, H.L., "Scaling Effects on Bubble Dynamics in a Tip Vortex Flow: Prediction of Cavitation Inception and Noise," Rep. 98007-lNSWC, 2000, Dynaflow, Inc., Jessup, MD. [1SJ Hsiao, C.-T., Chahine, G.L., Liu, H.L., "Scaling Effects on Prediction of Cavitation Inception in a Line Vortex Flow," to appear in ASME Journal of Fluids Engineering, 2002. [16] Hsiao, C.-T., Chahine, G.L., "Numerical Simulation of Bubble Dynamics in a Vortex Flow Using Navier-Stokes Computations and Moving Chimera Grid Scheme," 4t International Svmposium on Cavitation, Pasadena, CA, June 20-23, 2001. t17] Johnson, V.E., Hsieh, T., "The Influence of the Trajectories of Gas Nuclei on Cavitation Inception," Sixth Svmposium on Naval Hydrodynamics, 1966, pp. 163-179.
[18] Ligneul, P., Latorre, R., "Study on the Capture and Noise of Spherical Nuclei in the Presence of the Tip Vortex of Hydrofoils and Propellers," Acustica. Vol. 68, 1989, pp. 1-14. El91 Ligneul, P., Latorre, R., "Study of Nuclei Distribution and Vortex Diffustion Influence on Nuclei Capture by a Tip Vortex and Nuclei Capture Noise," ASME Journal of Fluid Engineering, Vol. 115, 1993, pp. 50~507. [20] Plesset, M. S., "Dynamics of Cavitation Bubbles," Journal of Applied Mechanics, Vol. 16, 1948, pp. 228-23 1. t21] Roe, P. L., "Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes,: Journal of Computational Physics, Vol. 43, 1981, pp. 357-372. t22] Shen, Y.T., Jessup, S.D., Cowing, S., "Tip Vortex Inception Logarithmic Scaling," NSWCCD-50-TR-2001/060, October 2001a, Naval Surface Warfare Center, Carderock Division, West Bethesda, MD. [231 Shen, Y.T., Chahine, G.L., Hsiao, C.-T., Jessup, S.D., "Effects of Model Size and Free Stream Nuclei on Tip Vortex Cavitation Inception Scaling," 4th International Symposium on Cavitation, Pasadena, CA, June 2~23, 2001b. [24] Taylor, L.K., Pankajakshan, R., Jiang, M., Sheng, C., Briley, W.R., Whitfield. D.L., Davoudzadeh, F., Boger, D.A., Gibeling, H.J., Gorski, J., Haussling, H., Coleman, R., and Buley G., "Large-Scale Simulations for Maneuvering Submarines and Propulsors,", AIAA Paper 98-2930, 1998. [25]van Leer, B., "Towards the Ultimate Conservative Difference Scheme. V. A Second Order Sequel to Godunov's Method," Journal of Computational Physics, Vol. 32, 1979, pp. 101- 136. [26jVanden, K., Whitfield, D. L., "Direct and Iterative Algorithms for the Three-Dimensional Euler Equations," AIAA-93-3378, 1993. [27] Wilson, R., Paterson, E., and Stern, F. "Unsteady RANS CFD Method for Naval Combatant in Waves", 22n~ Svmposium on Naval Hvdrodvnamics, Washington, D.C., 1998, pp. 183-197.
DISCUSSION S. Cordier Bassin d'Essais des Carenes, France The results shown of a novel modelization of bubble dynamics are very interesting The importance of nuclei size and density in cavitation (surface and vortex) is well put in evidence. This effect has been demonstrated in experiments in the GTH and reported in the literature (Fruman and Briancon). Do you have any plans to validate your results against this type of data? AUTHORS' REPLY We are well aware of the work you mentioned. The results shown in our paper are focused on single bubble dynamics. Validation of our current scheme requires a well-controlled experiment in which the initial nucleus conditions are known and the nucleus is followed in its motion. To study the effect of nuclei size and density on the scaling problem, a statistical type model which allows one to randomly distribute the nuclei according to a known nuclei size distribution is required. Since this will require a large number of simulations with a significant nuclei population, the SAP spherical model, described in the paper, seems to be a good choice for such a study that we are presently conducting. We will be presenting such results in the future and will try to validate our results with the experimental data you suggested assuming you could provide us with the nuclei size distribution in those tests. DISCUSSION R. Gornstein Navatek Ship Ltd., USA Can your method be extended to surface flows (sheet cavitation) as well as vortex flows? AUTHORS' REPLY Actually, our numerical scheme was initially developed to solve the sheet cavitation problem (Chahine and Hsiao 2000~. We then extended it by incorporating a moving Chimera grid scheme to enable the simulation of the dynamics of moving bubbles in a tip vortex flow. The same scheme, of course, can be used to further improve the simulation of sheet cavitation. However, the current scheme is limited when the interface becomes multi-connected. We are working on further improving this aspect.
