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.
Finite Difference Analysis of Unsteady Cavitation on a Two-Dimensional Hydrofoil A. Kubota, H. Kato and H. Yamaguchi University of Tokyo Tokyo, Japan Abstract The authors present a new cavity model, which is named a Bubble Two-phase Flow (BTF) model, to explain the interaction between viscous fluid and bubble dynamics. This BTF cavity model treats the inside and outside of a cavity as one continuum by regarding the cavity as a compressible viscous fluid whose density varies widely. Navier- Stokes equations including cavitation bubble clusters are solved in finite difference form by a time-marching scheme. The growth and collapse of a bubble cluster is given by a modified Rayleigh's equation. Computation was made on a two-dimensional flow field around a hydrofoil NACA0015 at an angle of attack of 8.0 deg. The Reynolds number was 3x105. The computational results showed the occurrence of the attached cavity at the foil leading edge. Furthermore, the present results predicted the generation of cavitation clouds with large-scale vortices. The newly proposed BTF cavity model is very flexible and promising. 1. Introduction When the pressure of a liquid is reduced at constant temperature by either static or dynamic means, a state is reached at which vapor or vapor- filled bubbles exist. This phenomenon is called cavitation [1]. Cavitation results in performance deterioration of hydraulic machinery, generation of noise, vibration and erosion. This is because cavitation is a dynamic phenomenon, as it concerns itself with the growth and collapse of cavities. In spite of many excellent studies, the actual structure of cavitation is not 667 yet fully understood. In order to accurately predict the generation of noise or onset of erosion, it is particularly important to illuminate the unsteady structure of cavitation. Vortex cavitation is often observed in the flow downstream of attached cavitation. It is caused by vorticity shed into the flow field just downstream of the cavity. Such vortex cavitation generates a large cavitation cloud under certain conditions. The vortex cavitation impinges on the body and its subsequent collapse results in erosion [2]. In previous work, the authors performed an experimental investigation of the unsteady structure (velocity distribution) of cloud cavitation on a stationary two-dimensional hydrofoil using a conditional sampling technique [3]. It was found that the cloud cavitation observed in the experiment was a large-scale vortex with many small cavitation bubbles. Consequently, the importance of the interaction between large-scale coherent vortices in the flow field and cavitation bubbles was recognized. Much theoretical work also has been done in order to obtain a better understanding of the physics of cavitation phenomenon. Researchers have continually developed new models of cavitating flow based mainly on the assumption that the flow is irrotational (inviscid). For example, Tulin [4][5] proposed small-perturbation (linearized) theory for the case of a supercavitating hydrofoil. He has treated the cavity as a single vapor film and assumed that the pressure inside the cavity is constant. This is a sort of macroscopic analysis of cavitation. Using this model, many researchers have gradually improved the
calculation methods [6][7][8][9][10][11]. The single vapor film model is now well established. It can also predict macroscopic cavity characteristics fairly well. However, they all dealt with only steady cavitation. Furness and Hutton [12] treated the case of an unsteady attached cavity on a stationary two-dimensional body by the singularity method. This calculation result showed unsteadiness of the cavity surface and a reentrant jet [1]. Their methods, however, could not predict the behavior of a detached cavity after the attached cavitation splits into two parts. Tulin and Hsu [13] and van Houten [14] have solved the unsteady cavity problem on a periodically oscillating hydrofoil. Their method also could not predict the generation of detached cavitation clouds. This is because of a limitation of the cavity model, which treats the cavity as a single vapor film where pressure is constant. Therefore, a new model of cavitation is required to theoretically study the breakoffs of attached cavitation and cavitation clouds. How to model the cavity trailing edge, where the cavity collapses, is the most difficult problem for the above- mentioned single film, constant pressure cavity model. When one observes actual cavitation, it is found that the sheet cavity splits into minute bubbles with vortices in the end region. Then the bubbles collapse. One can often observe many vortex cavities in this region even if the sheet cavity is stable. Van Wijngaaden [15][16], March [17], Chahine and Lie [18], d'Agostino and Brennen [19] and others have studied the dynamics of bubble clusters. The bubble cluster is a kind of microscopic modeling of cavitation. However, these studies treated only a cluster of collapsing bubbles under given conditions. Hence, they could not answer how the unsteady attached cavity sheds cavitation clouds. Highly vertical fluid motion such as a cavitation cloud is often observed downstream of a cavity. Experimental observation shows a close relationship between large-scale coherent vortex and cavitation [3]. Therefore, it is necessary to develop a theoretical elucidation of the mechanism that generates the large-scale vortex structure. There is a strong interaction between the large-scale vortex and cavitation. The occurrence of attached cavitation yields boundary layer separation. The separated shear layer rolls up, thus turning into a large- scale vortex [20]. On the contrary, the large-scale vortex yields a low pressure region at its center. In the low pressure region, bubbles grow and remain. Existing cavity models are powerless to explain the nonlinear vortex dynamics of the above problem since they assume inviscid flow. In this paper, the authors will propose a new cavity model that can explain the interactions between vortices and bubbles. The authors call this new model a Bubble Two-phase Flow (BTF) cavity model. In a macroscopic (coarse-grained) view, this model treats the cavity flow field phenomenologically as a compressible viscous fluid whose density varies greatly. It can treat the inside and outside of the cavitation as a single continuum and hence it can express the detached cavitation clouds. In a microscopic view, this model treats cavitation structurally as bubble clusters. By coupling these two views, the BTF cavity model can clarify the nonlinear interaction between macroscopic vortex motion and microscopic bubble dynamics. The authors have developed a program code SACT-III (Solution Algorithm for Cavitation and Turbulence, version III) to solve the BTF model equations using the finite difference method. They will apply this program to cavity flow around a stationary two-dimensional hydrofoil in order to verify its ability. 2. Formulation An attached cavity, which is formed at the leading edge of a cavitating hydrofoil, collapses at the mid-portion of the hydrofoil. It is also well known that the attached cavity oscillates cyclically within a certain range of cavitation number. This unsteady cavity sheds a cavitation cloud in each cycle [3]. The front part of the attached cavity is a film of vapor where pressure is constant. At its rear part, the vapor film splits up into tiny bubbles. A large-scale vortex caused by the cavity rolls up the bubbles, thus generating a cavitation cloud. There are two cavity types in a microscopic view. As shown in Figure 1, these are the vapor film and the bubble cluster (with vortices). 668
/ / . ;.; ~.\ CAVITY SURFACE // N[~ ( ( ~ ° ~ o of ~ ) ~ At_ - I;igure 1 Two microscopic views of cavitation, vapor f i lm types (above) and bubble cluster types (below). There is, however, a third type of cavitation. This is traveling cavitation or bubble cavitation [1]. This phenomenon is not treated in this paper. 2.1 Macroscopic Modeling In the macroscopic view, the Bubble Two-phase Flow (BTF) model treats the inside and outside of the cavitation as one continuum. This is because it regards the cavity flow field as a compressible viscous fluid whose density varies greatly. According to this phenomenological modeling, contour lines of void fraction (volume fraction of cavities) express the shape of the cavity as shown in Figure 2. Governing equations of the macroscopic flow field are as shown below. The equation of continuity is 50e+ Ad pv)=o , (~) p=(l-f g) AL 669 ~ I CROSCOP I C V I EW\ At\ 1~ , \ LOCAL HOMOGENEOUS MODEL / \ (LHM) ~ CAVITATION CLOUD _ MACROSCOP I C V I EW COMPRESSIBLE VISCOUS FLOW l~igure 2 Modelling concept of Bubble Two-phase Flow (BTF) cavi ty model. where t, pv, v and p are time, mass flux vector(pu, pv, pw), velocity vector (u,v,w) and density of the mixture, respectively. The conservation equation for momentum: Navier-Stokes equation is 8(G~+V(pvv)=-VP+Rlep{v2v~lv(v·v)} ,(2) where P is the pressure in the mixture, ~ is the viscosity of the mixture, and Re is the Reynolds number. This equation is in conservation form [21]. The nondimensionalized quantities based on the uniform flow velocity and a reference length have been employed in the above two equations. The BTF cavity model assumes that a fluid of variable density replaces the water-vapor mixture. The density of the water containing bubbles (mixture) is defined as follows: (3)
where PL is the water density and fg is the local void fraction. The mass and momentum of the vapor are ignored, since they are small compared with those of the liquid. The actual ratio comparing the density of vapor to that of liquid is of the order 10-4. The change of liquid mass due to the phase change is also ignored. The viscosity of the mixture is assumed to be as follows: P=(1-fg)PL+fgI`G ' (4) where ILL is the water viscosity and FIG is the vapor viscosity. 2.2 Microsconin Mod ~ 1 i n ~ - Local Homogeneous Model - To calculate the macroscopic flow field, it is necessary to know the local void fraction function fg(t,x,y~z). The greatest problem is to develop a model that gives the relationship between the flow field condition and the void fraction. The present BTF model treats cavitation microscopically as bubble clusters. This is because bubbles play important roles in the cavity inception [22]. Furthermore, one of the main purposes of the BTF model is to study the mechanism of the cavity collapse. The sheet-type cavity splits up into tiny bubbles there. This structural microscopic model cannot be basically applied to the vapor film type cavitation. However, the type of the microscopic model has less effect on the macroscopic model when the void fraction is high. The present BTF model introduces a Local Homogeneous Model (Lam), which is a sort of Mean Field Approximation (MFA), for simplicity. It treats the cavity as a local homogeneous cluster of spherical bubbles as shown in Figure 2. Bubble number density and a typical radius are assumed locally. This typical bubble radius is obtained from the growth-collapse equation of the bubble cluster. This equation is derived from the growth-collapse equation of one spherical bubble (Rayleigh's equation). The LHM gives the local void fraction fg by coupling the bubble density and the typical bubble radius as follows: ~ =n.4~R3 , where n is the bubble number density and R is the typical bubble radius (see Figure 3). In this paper, the bubble density is assumed to be constant all over the computational domain though real cavity flows have a distributed bubble density. This is because it is too difficult to formulate coalescence and fragmentation of the bubbles. , ~ if. _ ~ I.. ..~... .N ... . ..., . .. .. ....... ..... ...... ... ,) age:-.. \ //, ''.: f = n 4lrR3 p = ( l~fg ) PI, Figure 3 Local homogeneous microscopic cavity model (LHM). The LHM assumes that cavities form in the shape of spherical bubbles. The bubbles remain separate and distant enough from each other so that their shapes remain spherical. Interaction between bubbles occurs through the local pressure that develops in the liquid as bubbles grow. Lord Rayleigh originally derived the equation of radial motion (growth and collapse) of an isolated spherical bubble in a homogeneous infinite medium [23]. This equation is widely known as Rayleigh's equation. It takes the following form, neglecting the effect of surface tension and viscous damping, Rd2R+3(dR)2=Pv P (6) where Pv is the vapor pressure. In this study, the vapor pressure is assumed constant. This is because the behavior of the bubble is nearly isothermal and gas inside the bubble is also ignored. The finite-difference method was employed in SACT-III. In this method, a continuous domain is discretized into finite grid points. Hence, an (5) interaction between individual bubbles within the grid spacing must be considered to eliminate the effect of the computational grid. This effect is a sub-grid-scale (SOS) bubble interaction. 670
Next, let us consider the SGS bubble interaction of the LHM, deriving analytically the equation of motion of the bubble cluster. ~ MA r ~~ r Figure 4 Sub-grid-scale (SOS) bubble interaction mode 1. As shown in Figure 4, we consider the influence of the other bubbles which exist inside of the distance Or (=grid spacing). The total velocity potential due to the other bubbles at the origin O is i 1 dri 2 >: ri dtRi (7) when Ri<<ri. From the local homogeneous assumption, i dr. V(~ r1 dtlRi2)=0, and Ri=R . (8) (9) The following equation is therefore obtained by adding the time derivative of Equation(7) to the original Equation(6): l 1 dR 2 d2R 3 dR 2= v ·(10) to ri dtR )+Rdt2~2(dt) p The number of bubbles which exist inside the sphere of radius or is nuder . (11) = Bt(dRR2n.2~r2) =2~r2(nR2d R2+dtR2dR+2nRdR) .(12) Combining Equations(10) and (12) and replacing the time derivative with Ot {'~(v~v)} , we obtain the LHM's equation of motion as follows: (1+2~r2nR)RD 2~(3+4~Ar2nR)(DDi<)2 +2~r2DnR2DR= v Referring to Equation(13), it is found that the effect of the other bubbles decreases with a decrease in fir. This is a preferable characteristic in solving the present problem. This is because the other bubbles do not affect the referred bubble if the grid interval was zero. Quantities in the above equations have been nondimensionalized based on the uniform flow velocity Up* and a reference length d*. Hence: P =~L UP 2p, t =td*, v =U~ v, c<> * * * * ~ =PL I' ~ AL I' n = no, R =d R. d * * * * * * x =d x, y =d y, z =d z, where * denotes dimensional values. In the following computation, the chord length of a hydrofoil has been chosen as d*. The Reynolds number Re, pressure coefficient Cp and cavitation number are defined as follows: U d U d ~ Re= ~ = ~ L VL AL P -P P 2PL U Pm PV 1 * *2 -2rL Us Then the first term of the left-hand 3. SACT-III Program side in Equation(10) becomes t(~: 1 dRR2)=~3t(dRR2~: rl) t(dRtR2n~ 14~r2dr) 671 (15) The program SACT-III is the third version of the SACT series. The purpose of the SACT series is to study theoretically the unsteady structure of cavitation using the BTF model. Over several years, we developed a program SACT-II (SACT, version II; two-
dimensional rectangular cell version). The results computed by SACT-II explained cavity formation caused by large-scale coherent vortices behind a rectangular obstacle on a wall [24]. The present SACT-III employs the finite difference method in the body fitted coordinates to solve the governing partial differential equations given in the preceding section. This section explains the computational procedure and the finite difference scheme of SACT-III. The computational procedure is basically the same as the Marker-and-Cell(MAC) method [25] except for the use of a regular mesh system, instead of the staggerd mesh system. 3.1 Quasi-Poisson Equation for Pressure By taking the divergence of the Navier-Stokes equation(2), the following Poisson equation for pressure is given: where ~(~v,v)=-V[V(pvv)+Ilep{V2v+1V(V v)}] .(17) Substituting the continuity equation(1) into Equation(17), 2 2 V P= ~ 2p+~(pv,v) . ( at From Equations(3) and (5), which represent the basic assumptions of the LEM, v2p=_V(~ )+~(~v,v) =-~t{V(~)}~(~v,v) , (16) 18) ~=(l-n.3~R3)~L . (19) To simplify the LHM, n=constant is assumed in this study as mentioned in section 2.2. By differentiating Equation(19) twice with respect to t, we have the following equation: ^=-~L4n~{R a 2+2R(~t) } . (20) and From the LHM's equation of motion(13), ~ R=s7(pv,v,~R,R)+~(P) , (21) where ( pv,v, aR,R) =-{2(v.V)3R+(v.V)(v.V)R} _+4~Ar2nR BR 2 (22) 2 2 { t+(v.V)R} (1+2~^r nR)R ~ _ 2~Ar2R2 {8n+(v~v)n}{8nt+(v~~)R} , (l+2~r2nR)R ~t ~ 672 P. -P *(P)= 2 (23) (1+2~Ar nR)RPL If n=constant, the last term of the right-hand side becomes zero in Equation(22). Substituting Equations(20) and (21) into (18), we obtain the quasi-Poisson equation for pressure including the motion of the bubble cluster as follows: V2P+~'(P)=~'(pv,v,~R,R)+~(pv,v) ,(24) where ~'(P)=-pL4n~R ·~(P) , (25) a~ (pv,v,~3R,R) PL4n7tR { g! ( PV, V, a8R, R ) + 2 ( BR ) 2 } 26) The left-hand side in Equation(24) is approximated by the second-order finite differencing scheme. As a consequence, simultaneous equations of pressure P are obtained if the right-hand side is given in Equation(24). SACT-III solves the simultaneous equations derived from Equation(24) with a point successive relaxation method. Equation(24) is equivalent to the normal MAC method's Poisson equation of incompressible flow [25] when ~' =97' =0 They are always set to zero for non- cavitating conditions. If the mixture is filled with liquid, Equation(13) cannot be solved since the bubble radius R becomes zero. If the mixture is filled with vapor, Equation(2) cannot be solved since density of the mixture becomes zero. Hence, when the void fraction fg isless than fgmin(>0 °) or more than f~max (<1.0), the bubble radius becomes f~xed. Then, ~ and a'' are also set to zero. In the following computations, fgmin and fgmax are set as follows: i =n 4~Ro3 ' f =0.95 , where Ro is the initial bubble radius. 3.2 Numerical Methods Equations(2) and (13) are time- integrated with and Euler explicit scheme using the value of pressure P obtained by solving Equation(24). To solve a high-Reynolds-number flow, each nonlinear term in Equation(2), for example ~3x(~uv), was approximated with
the fourth-order centered finite- differencing scheme with the fourth- derivative term: 4 ul. a4(pv)^x3 . (27) ax The fourth-derivative term plays an important role in stabilizing the calculation. Physically, fourth- derivative term means shorter-range diffusion compared with the second- derivative viscous term [26]. The fourth-derivative term consequently stabilizes the computation without decreasing the Reynolds number and introducing any turbulent models. The universal availability of this fourth- Downst ream derivative term nas not been ascertained yet for the turbulent flow calculation. Boundary However, SACT-III introduces no turbulent model since there exists no established one for the two-phase (cavity) flow at present. All the other space differential terms in Equations(2) and (13) are approximated with the second-order centered differencing scheme. Equation(13), however, has no spatial diffusive term for bubble radius R and its time-derivative. The second- derivative term is accordingly added in equation (13) to eliminate the instability of the nonlinear terms. For example, lu1 ·B Rex ax is added to u8R . This term means diffusion. To compute the high-Reynolds-number flow around a body of arbitrary shape, it is convenient to use body-fitted coordinates through coordinate transformation. Figure 5 shows the grid system for the present problem of flow around a two-dimensional hydrofoil. This system is called a C-type grid [27]. The connected physical (x,y,z) domain around the hydrofoil is mapped onto the rectangular computational (\,q,;) domain. Here the pair of planes forming the branch cut are both on the same plane of the transformed region. The surface of the body is also mapped on the same plane with the branch cut. A regular mesh system is employed. Velocities, pressure and bubble radius are given on the grid points. As shown in Figure 5, the uniform flow boundary PLOW ~ ~ = PHYS I CAL DOMA I N Outer-flow Boundary Side Boundary Downstream Boundary ... . j ~ Wall Boundary ~ L ~ _ _ _ ~ Branch Cut TRANSFORMED DOMA I N Fi Sure 5 C-grid system around a two-dimensional hydrofoi 1. conditions are imposed at the outer-flow boundary. These are: u=l, v=w=O, Bu_8v_8w_O B~ B9~ all- , R=R , ~R=0 P=0 . O at ' (29) 28) At the downstream and side boundaries, the boundary conditions of the zero-th order (zero-gradient) extrapolation are imposed. At the branch cut boundaries, the periodic boundary conditions are imposed. At the wall boundary, the following boundary conditions are imposed. u=v=w= 0, BP=BR=0 (30) First order (linear) extrapolation in (t,Q,4) domain is used for the velocities to calculate the nonlinear terms of Equation(2). The terms are approximated by the fourth-order centered differencing scheme with the fourth-derivative term. 4. Computational Results and Discussion 4.1 Condition of Computation 673 A hydrofoil section with a simple mathematical configuration, NACA0015 [28], was chosen for the computation.
The computation was performed at angles of attack ~ of 0.0, 8.0 and 20.0 degrees [29]. However, this paper discusses only the results at ~ =0.0 and 8.0 degrees. The Reynolds number Re was 3x105 in all the computations. It is based on the uniform flow velocity and chord length of the hydrofoil The viscosity ratio G/pL was 9.12x10-3. The computation at ~ =0.0 deg. was performed only for non- cavitating conditions to evaluate numerical accuracy. Experimental cP observation at ~ =8.0 deg. shows laminar separation without bursting near the leading edge. For cavitating conditions, an attached type cavity accordingly occurs near the foil leading edge. The hydrofoil was accelerated from u=0 to the steady speed of 1 for T=O-1. The cavitation number was also decreased gradually for T=2-3 so as to compute stably for cavitating conditions. The time increment t was determined in each computational step to keep the Courant number less than 0.25. The relaxation factor was set to be 0.8 when Equation(24) was solved with the successive relaxation method. However, it was reduced to 0.3 for cavitating conditions. The convergence condition of the pressure computation is as follows: max(lAPI) < 0.001 (for non-cavitating conditions) max(lAPI) < 0.01 (for cavitating conditions), where UP is the residue of pressure in an iterative calculation. Figure 6 C-grid system around NACA0015 hydrofoil at a=O. Odeg, lOl(~)x31~ ~ x3~. Preceeding the computations for cavitating conditions, it is necessary to evaluate numerical accuracy of SACT-III for non-cavitating conditions. The angle of attack was 0.0 degrees. Figure 6 shows the c-grid system. The grid was uniform in the spanwise section. The number of grid points was lOl(~)x31(~)x3(~). The distance between the trailing edge and the upper or lower boundaries was 1.2. The distance between 674 the trailing edge and the downstream boundary was 3.0. The shape of the front part of the outer-flow boundary was oval as shown in Figure 6. A direct numerical method [30] was used for the grid generation. The minimum grid spacings at the leading and trailing edges were 0.70xlO 4 and 1.54x10-4, respectively . The minimum spacing was about one twelfth of Re~0 5 at the trailing edge. -2.0 - _ -1.0- _ · Experiment (Re=6 x 10 5) 06 Present Calculation (Re=3x 105) Hess-Smi th Method Hi th Boundary Layer Calculation (Re=3 x lOs) ~o57 1 Am_ I', , , 1 `49E ,^~ ,______ ~ Figure 7 Foil surface and wake pressure, NACl\OO 1 S; ~ =0. Odeg. ; T=2. O. Figure 7 shows the pressure distributions on the foil surface at T=2. The circles show the distribution of Cp on the foil surface. The triangles show the distribution of Cp in the wake, i.e., on the branch cut shown in Figure 6. The computational results are compared with the computation using the Hess-Smith method [31] and the measurement at Re=6xlO5 by Izumida [32]. The Hess-Smith method is a sort of numerical solution method of potential flow based on the boundary element method. The foil shape was modified by adding the computed displacement thickness of the foil surface boundary layer [33]. The laminar boundary layer was computed using Thwaite's method, with the empirical constants derived by Curle and Skan [34]. The length of the separation bubble was 150 times as long as the momentum thickness at the laminar separation point. The Head's entrainment method modified by Cebeci [35] predicted the turbulent boundary layer development. The turbulent separation was predicted to occur when the form factor H12 exceeds 2.1. The computed pressure distribution agrees very well
with the others as shown in Figure 7. Furthermore, the pressure coefficient is 0.982, which is almost equal to 1.0, at the front stagnation point. As a consequence, the present numerical method has good accuracy when the grid system is fine enough. The computed pressure distributions disagree with the experimental result when the coarser grid systems are used [29]. Convective Term _~ ~ ~ DMAX=14.4 Fourth-derivative Numerical Dissipation Term Pressure Term Truncation Error of Pressure Term Figure 8 Distribution of convective, fourth- derivative, pressure and its truncation O error terms in x-momentum equation, AL. NACA0015; a=O. Odeg. ;Re=3. Ox 105;T=2. O. --- Figure 8 shows the distributions of the convective term, the fourth- derivative term, the pressure terms and its truncation error in the x-momentum equation. The contour interval is 1.0. DMAX shows the maximum. The maximum of the fourth-derivative term is not negligible compared with the convective and pressure terms. However, it is distributed only near the foil surface around the leading and trailing edges. This result shows that the effect of the fourth-derivative term is local but important to stabilize the computation. The truncation error of the pressure term is sufficiently small compared with the main differences terms. 4.2 Unsteady Attached Cavity Figure 9 shows a close-up of the grid system around the NACA0015 hydrofoil. The angle of attack was 8.0 deg. The number of grid points was lOlx31x3. It was the same as that at ~ =0.0 deg. The minimum spacing was 0.70x10-4 at the Figure 9 Close-up of the grid system around NACA0015 hydrofoil, a=8.0deg., 101~) X31~7~) x3~. trailing edge. It was almost the same as that of the grid system at ~ =0.0 deg. Figure 10 shows the time-averaged velocity vectors from T=2 to 4 for non- cavitating conditions. The computation was performed stably. The boundary layer separates at X=0.74 on the back side. Instantaneous velocity vectors, however, show unsteady vortex shedding from the foil trailing edge region. In the other region, the flow is almost steady. No separation occurs near the leading edge. SeparationPoint (X=0.74) ~_~ _~` _ 675 Figure 10 Time-averaged (T=2~4) velocity vectors around NACA0015 hydrofoil, a=8.0deg.; Re=3. 0 x 10 s. Figure 11 shows the chordwise distribution of boundary layer displacement thickness on the foil surface. The present result agrees well with the boundary layer calculation at the front part of the back side. However, it cannot predicte a laminar separation bubble because of the insufficiency of the grid points in the }-coordinate. The separation point is further upstream than in the boundary layer calculation. On the contrary, the agreement of the separation point is good on the face side. However, the boundary layer is thicker near the leading edge. Figure 12 shows an example of the velocity profiles using wall variables. The present profile closely follows the low of the wall. The present method can consequently express the turbulent boundary layer without any turbulent models.
~ 1 a - ut 20: 10 o BACK SIDE Present Cal. -- Boundary Layer Ca 1. ·: Separat ion Po i nt t: Separ at i on Bubb l e > T~rh~ 1 An t ~ 7 FACE SIDE However, the shape of the pressure distribution agrees well with the experiments. Hence, the disagreement of pressure hardly affects the nature of the unsteady cavitation. ////// ~ Pv- P p x _ = -O. 1 GO\ L 1` ~ A: non ~ r3 = to B: n3'T ~ r3 = 1 1 C: n-n ~ r3 = 10 *1 C ~ |B | D: n3'T ~ r3 = 100 Present Cal. - - - Boundary Layer Ca 1. / P v- P T i ~ e t p ·: Separation Point ~ Gin , x 0 0.5 +0. 1 3.2 |.0 t.8 _ 2.' Figure 11 Boundary layer displacement thickness .= 2.2 distribution on NAChOOlS hydrofoil. ~ ct=8.0deg. ;2e=3.0X105. ~ 2.0 _ I.~ ~ l.. nut ye '/Z/571ny++5.5 . :Present Cal. (I=57) Of = 5. 91x 10-3 v*O= 8. 07x 10-2 Figure 12 Velocity profile in boundary layer on NACA0015 hydrofoil using wall variables, ~ =8.0deg.;Re=3.0X 105;I=57(X=0.0827~. While the velocity profiles are predicted well, the pressure is not predicted quite so well. The present lift coefficient CL is only about 58% that obtained by the experiment. ~ A 4 ~~ 3 3 B : n 3 7r ~ r 3 = 1 - C: non ~ r: ~ D 4 a 3 - = 100 0 1 2 3 ~ 5 e 7 8 ~ 10 Tise t Figure 13 Influence of the SGS bubble interaction model on bubble collapse (above) and growth (below). Before analyzing the cavitating flow field, we should check the effect of the SGS bubble interaction model. Figure 13 shows the effect of the SGS bubble interaction when the ambient pressure P has changed stepwise with the amplitude of +(Pv~P) /PL=0 .1. n(4/3)= ~r3 is the number of bubbles within the grid scale a r. The convective terms were neglected in Equation(13). The Runge-Kutta-Gill method was used to solve the of differential Equation(13). The bubble radius is nondimensionalized using its 676
initial value. The case A is the calculation of isolated bubbles, i.e., with no-bubble interaction effects. As shown in this figure, presence of other bubbles delays growth and collapse of the bubbles. Fujikawa et al. [36] carried out a theoretical analysis on the interaction between two bubbles. Their result shows that the collapse is delayed when the two bubbles have the same radius. The present result shows the same tendency. In the following computation, the grid scale of the SGS bubble interaction model is assumed as follows: r=(~)2 where g' =X.,y~-x~y em (31) (two-dimensional Jacobian). This is because the flow structure is two-dimensional. Accordingly, the SGS bubble interaction effect is independent of the grid spacing in the spanwise direction. /; Vertical Exaggeration of 4:1 Figure 14 Void fraction contours around NACAOOlS hydrofo i l, cr = 1. 5; ~ =8. Odeg. ; Re=3. 0 x 10 5 - the contour interval is 0.1 except for the outermost line. Figure 14 shows contour lines of void fraction at ~ =1.5. The initial bubble radius Ro and the bubble number density n are 4x10-4 and lx1O6, respectively. In this paper, we define a cavity as a region where void fraction is more than 0.1. The bold contour lines are those of fg=O.1 in this figure. Contour lines are drawn at intervals of O.1 except for the most outer line of fg=O.01. As seen in this figure, thin steady attached cavity incepts smoothly. Figure 15 shows the time series contour lines of void fraction at 0.1 0.01 W\_ ~F~ \ T = ~ OF ~ I; igure IS Void fraction contours around NACA0015 hydrofoi l, ~ =1. 2; a =8. Odeg. ; Re=3. 0 x 10 5; the contour interval is 0.1 except for the outermost line. ~ =1.2. The void fraction is more than O.9 at the center of the cavity. The rear portion of the cavity oscillates cyclically. Not only does the cavity length change but the cavity itself rises up at its rear part. The unsteady characteristics computed here agree with the experimental observations of sheet- type cavitation [32]. It is therefore concluded that the BTF cavity model can express the features of sheet-type cavitation beyond its microscopic model, which is essentially suitable to the bubble cluster flow. Figure 16 shows the time-averaged pressure distribution on the foil 677
T I ME-AVERAGED CP-DISTRIBUTION CP (T=3. 0~ 5. 0) NIiCt0015; a =8O, ~ =1. 2 O . _1 n n 1.0 CL = 0. 5472 CD = 0. 0693 ........ it_ ~ L. Ed H.C -A- ~ - ' ~ 0.0 ~_0. · P.~v i t v A rea ~ - ~ ~` ,J,,,_~ ~ Figure 16 Time-averaged pressure distribution and void fraction contours around NACA0015 I; hydrofoil, a=1.2;a=8.0deg. ;Re=3.0X105, the contour interval is 0.1 except for the outermost line. surface and void fraction contour lines. The foil surface pressure where cavity exists is almost constant and equal to the vapor pressure (Cp=-1.2). However, a small pressure peak exists at the front part of the cavity. The pressure distribution is similar to the calculated result by a nonlinear free- stream line theory [10]. Furthermore, the time-averaged cavity shape is similar to the experimental observation of sheet-type cavity. These facts also suggest that the present BTF model can be applied to sheet type cavitation beyond its micro structure limitation. Vertical Exaggeration of 4:1 Contour Interval of 0.1 ~\ _ _ _ o °L i gu re 17 T i me-averaged pres sure coef f i c i ent contours and velocity vectors around NI\CAOOlS hydrofoi l, a =1. 2; a =8. Odeg. ;Re=3. 0 x 105, the contour interval is 0.1, the bold broken lines are void fraction contours of 0.1. Y r o.l r o.lo ~ I MU = U O t.0 0 1.0 OISPLRCEt1ENT Tl~lct~rlEss = 0.0288q DlSPLPCEMENr THICKNESS = 0.02112 hlOMft`TUb TFilCht~fSS ~ O.00q97 MOMfNTUM THICKNESS ~ 0.00557 FOnM FnCTORll112) = S.8O FORM FRCrOlltH12} = 3.79 o. 10 ~ r O 1=67 1 0.10 1=69 X=O. 4812 X=O. 5832 1 o 1.0 o 1.0 Figure 17 shows the time-averaged pressure contours and flow velocity vectors. The contour line of fg=O.1 (the bold broken line) agrees approximately with that of Cp=-1.2. This is the reason why the cavity is defined as a region where void fraction is more than 0.1. The reverse flow is observed clearly in a=1.2 the neighborhood of the end of the ~ - No-cavitation Bold lines in Figure 18 show time- averaged velocity profiles in the cavity wake boundary layer along the 1-coordinate. Fine ones are profiles DISPtnCEh1Et1T li1ICKt1ESS = 0.01~! OlSPLnCEMENT TFfICKNE'SS = 0.01507 MOl1ENTUM Tt1ICK'lESS · 0.00771 MOMENTUM THICKNESS ~ 0.00813 FORt1 fnCTOnt'll2, = t.f39 FORH F0CTOR{~12) = I.f3S I;igure 18 Time-averaged velocity profiles in cavity wake region, a =1. 2; a =8. Odeg. ;Re=3. 0 x 105. fine lines show the boundary layer velocity profiles for non-cavitating conditions. 678
~ 1 r=1.2 O. 03 O . O O. O to No-cavitation //// / / / I /1 / ,{ l Figure 19 Comparison of boundary layer displacement thickness distribution on the back side on 4.3 Cloud Cavitation Figure 20 shows void fraction contour lines at ~ =1.0 with a time increment of 0.2. The initial bubble radius Ro and the bubble number density n are 1x10-3 and 1x106, respectively. For this condition, unsteady cavities continuously grow and collapse. The highly distorted attached cavity sheds cavitation clouds cyclically (T=5.9 and 7.1), which soon collapse. This phenomenon agrees well with many experimental observations [32][39][40][41][42][43]. Figure 21 shows an example of the photographs of cloud cavitation [29]. The position of the cavity break off point agrees well with the computational result. Figure 22 shows velocity vectors around the foil. Overlaid bold broken lines are void fraction contours of 0.1. As mentioned before, they show instantaneous cavity shapes. As shown in this figure, the unsteady attached cavity sheds not only cavitation clouds but also vortices (see marks A, B and C). The experimental result has confirmed such vortex shedding phenomena [3]. The position of cavitation clouds, however, does not agree very well with That of shedding vortices. There are two ~ - - - - - - - 5 possible explanations for this NACA0015 hydrofoil. a =8.0deg.;Re=3.0X10 . for non-cavitating conditions. As shown in the section of 1=61 (X=0.2115), the flow inside the cavity is quite slow except near the foil surface. At the section of 1=64 (X=0.3375), strong reverse flow occurs. At 1=67 and 69 (X=0.4812 and 0.5832), the flow reattaches and inflection points exist in the velocity profiles. The measurement result downstream of a stable sheet cavity shows similar inflection points in the velocity profile [37][38]. As shown in this figure, the generation of cavity causes an increase in the boundary layer thickness behind it. Figure 19 shows the comparison of the chordwise distributions of the boundary layer displacement thickness for cavitating and non-cavitating conditions. For cavitating conditions, the displacement thickness decreases at the cavity collapsing region. It becomes minimum near the mid-chord, then it begins to increase again. This is also the same tendency that the experimental results show [38]. discrepancy. One is the phase delay of the cavity growth. The other is that the present cavity model does not consider nuclei convection towards the center of vortices. Figure 23 shows a close-up of the velocity vectors around the cavity. This figure elucidates the mechanism of cavitation cloud shedding. At T=5.5, a new separation vortex occurs at the cavity leading edge. Then it induces the flow toward the foil surface (T=5.7). Fluid density and pressure on the foil surface increase due to the impinging flow. It causes the cavity to break and tear off (T=5.9, separation of the cavitation cloud). The impinging flow turns into a jet along the foil surface. The jet sweeps away the cavitation cloud (T=6.1). This is the scenario of the generation of a cavitation cloud. 5. Concluding Remarks In this study, the authors presented a new modeling concept of cavitation called BTF (Bubble Two-phase Flow). In a 679
~- T = 5. T = 5. 5 \ it_ -I= ~~\ _~-U.01 - T = 6. 9 - CavitationCloud ~°~N \ T = 6. 1 ~ Figure 20 Void fraction contours around NACA0015 hydrofoil, a=1 0; a=8. Odeg. ;Re=3.0x 105, the contour interval is 0. 1 except for the outermost line. macroscopic view, this new cavity model treats the inside and outside of a cavity as one continuum. That is, it regards the cavitating flow field phenomenologically as a compressible viscous fluid whose density varies greatly. Contour lines of void fraction can express the cavity shape. In a microscopic view, a simple LHM (Local Homogeneous Model) is introduced. This is a kind of Mean Field Approximation. This structural microscopic model treats a cavity as a locally homogeneous bubble cluster. Assuming bubble density and a typical bubble radius, a local void fraction function is given. The BTF cavity model is significant in the following points: (l) The BTF cavity model can investigate the nonlinear interaction between large-scale vortices and cavitation bubbles, (2) The BTF cavity model can consider the effects of bubble nuclei on cavitation inception, 680 N: -,,,,! _ Hi_ ; Figure 21 Cavity appearance on NACA0015 hydrofoil, ~ =1. 30; a=8. Odeg. ;Re=3. Ox 105. (3) The BTF cavity model can express unsteady characteristics of cavitation. The BTF cavity model, therefore, includes three essential factors for cavitation. Those factors are pressure, nuclei and time. By examining the
- - - - - - ~: ~ - - - - - - ~ - - -. - - - - indispensable to improve the computation scheme of SACT-III to obtain qualitative agreement with experimental results. The accuracy improvement of numerical computations is one of the most important subjects at present. This research is partly supported by the Grant-in-Aid for Co-operative Research of the Ministry of Education, Science and Culture. All the computations in this study were carried out on the HITAC M-680H computer at the Computer Center of the University of Tokyo. The authors express cordial thanks to Mr. M. Maeda and Mr. M. Miyanaga for their help. They also express their appreciation to Prof. H. Miyata, Dr. Y. Kodama, Dr. N. Baba, Prof. H. Ohtsubo and Mr. J. Butler for their encouragement and support. References 1. Knapp,R.T. et al., Cavitation, McGraw-Hill (1970). 2. Hutton 9 S.P., Proc.Int.Symp.Cavitation, Vol.1, Sendai, Japan, pp.21-29 (1986). 3. Kubota,A.et al., J. Fluids Eng. Trans of ASME, Vol.111(6) (1989). 4. Tulin, M.P., Proc.NPL Symp. Cavil Hydro., Paper No.16, pp.1-19 (1955). 5. Tulin, M.P., J. Ship Research, Vol.7, No.3, pp.16-37(1964). 6. Wu, T.Y., J. Fluid Mech., Vol.13, Part 2, pp.161-181 (1962). 7. Nishiyama, T. et al.,Tech. Rep.Tohoku ,_,_,_~=,=,_, University, Vol.34, pp.123-139 (1969). _ 8. Furuya, O., J. Fluid Eng., Trans. ASME, Vol.71, pp.339-359 (1975). 9. Furuya, O., IAHR 10th Symp., Tokyo Fi gure 22 Veloci ty vectors around NACAOO15 hydrofoi 1, pp.221-241 (1980). a=l.O;a=8.0deg.;Re=~.OXl05,the bold brokenlo. Yamaguchi, H. et al.,Conf.on Cavil Iines are void fraction contours of 0.1. Edinburgh, IME, pp.167-174 (1983). 11. Lemonnier, H. et al., J. Fluid Mech., Vol.195, pp.557-580 (1988). 12. Furness, R.A. et al., J. Fluid Eng., Trans. ASME, pp.515-522 (1975). 13. Tulin T.P. et al., 13th Symp.on Naval Hydr., Tokyo, pp.107-131 (1980). 14. van Houten, R.J., 14th Symp.on Naval Hydr., Session V, pp.109-158 (1982). 15. van WiJngaaden,L.,Proc.llth Int.Cong Appl. Mech., Springer, pp.854-861 (1964) 16. van WiJngaaden, L., J. Fluid Mech., Vol.33, Part 3, pp.465-474 (1968). 17. M0rch, K.A.,Cavitation and Polyphase Flow Forum - 1981, pp.1-10 (1981). 18. Chahine, G.L. et al., J. Fluid Mech. Vol.156, pp.257-279 (1985). 19. d'Agostino, L. et al.,J, Fluid Mech. Vol.199, pp.155-176 (1989). computational and experimental results carefully, the BTF cavity model was proven to be useful to theoretically study cavitation characteristics. Complicated interactions, in particular, between large-scale vortices and bubble dynamics were clarified. However, further effort must be devoted to improvement of the microscopic cavity model for development of the BTF cavity model. For instance, the pressure gradient of the vortex cavitation cloud attracts bubbles towards its center. We must therefore take account of the convection of bubbles and the slip between bubbles and liquid to predict the behavior of cavitation clouds more accurately. Moreover, it is also 681
~l ~ ~ f ~' ol _l ~.0 \ L 1 ''' ~ 1.01} 1 ~ ~ , - ~ /, . ~' W\, \,,". ~'.,, , , , X/C 0 0.1 0.2 0.3 0.4 ~ / ~ , `\ _ '_~ ~ _ ~/ - _ \ 4-_ ~ ~ \1_ \. ~ _ ,. _ _ 74\\ _ /71 ~ _ _, ~Oo \ \ ,` .0t ,.90 ~ ~ I ' ' X/C 0 O.1 0.2 0.3 0.4 Vertical Exaggeration of 4:1 Figure 23 Close-ups of velocity vectors around NACAOO15 hydrofoil, a=1.0; a=8. Odeg. ;Re=3.0X lOs, the bold broken lines are void fraction contours of 0.1. 20. Kiya, M. et al., J.Fluid Mech., Vol.154, pp.463-491 (1985). 21. Roache,P.J., "Computational Fluid Dynamics", Hermosa Publishers(1976). 22. Kodama, Y. et al., J. Fluids Eng., Trans.of ASME, Vol.103(4),pp.557 (1981). 23. Lamb, H., Dover Pub., "Hydrodynamics", pp.112 (1932). 24. Kubota, A. et al.,Theoretical and Applied Mechanics, Vol.36, University of Tokyo Press, pp.93-100 (1988). 25. Harlow, F.H. et al., Phys. Fluids, Vol.8, No.12, pp.2182 (1965). 26. Kawamura, T. et al., AIAA paper 84-0341 (1984). 27. Thompson,J.F. et al.,"Numerical Grid Generation",Elsevier Science Pub.(1985). 28. Abbott, I.H. et al.,"Theory of Wing Sections", Dover Publications (1958). 29. Kubota, A., Doctral Dissertation, Department of Naval Architecture, University of Tokyo (1988). 30. Kodama, Y., J. Soc. Naval Architects of Japan, Vol.164, pp.1-8 (1988). 31. Hess, J.L. et al., Progress in Aero. Science, Pergamon Press, Vol.8 (1967). 32. Izumida, Y. et al, Proc. 10th IAHR Symp, Tokyo, pp.169-181 (1980). 33. Yamaguchi, H. et al., J. Soc. Naval Architects of Japan, Vol.164, pp.28-42 (1988). 34. Curle, N. et al., Aero. Quart., Vol.8, pp.257 (1957). 35. Cebeci, T. et al.,"Momentum Transfer in Boundary Layer",Hemisphere Co(1977). 36. Fujikawa, S.et al., Proc. Int. Symp. on Cavitation, Vol.1, Sendai, Japan, pp.55-60 (1986). 37. Yamaguchi, H. et al., LDV and Hot Wire/Film Anemometry, pp.29-38 (1985). 38. Kato, H. et al., Proc. Vol.2 18th ITTC, Kobe, pp.433-437 (1987). 39. Kermeen, R.W., Hydrodynamic Lab., CALTEC, Pasadena, California, Report No. 47-5 (1956). 40. Wade,R.B. et al.,J.Fluid Eng.Trans. ASME, Vol.88, No.1, pp.273-283 (1966). 41. Alexander, A.J., Conf. on Cavitation, Edinburgh, IME, pp.27-35 (1974). 42. Shen, Y.T. et al., 12th Symp. on Naval Hydro., Washington, pp.470-493 (1978). 43. Franc, J.P. et al., J. Fluid Mech., Vol.154, pp.63-90 (1985). 682
DISCUSSION by F. Stern I would like to congratulate the authors on a very ~ nteresting paper which appears to present a new approach for unsteady cavitation. The authors discuss the fact that cloud cavitation is often periodic. In my own work on uns teady cavi tat ion (Stern,F.,"Comparison of Computational and Experimental Unsteady Cavitation on a Pitching Foil", J. of Fluids Engineering, Vol.111, No.3, September 1989, pp.290-299), close correlation was shown between the experimental cloud-cavitation shedding frequency and the predicted cavity natural frequency. Do the authors' results provide an estimate for the cloud-cavitation shedding frequency and how does it compare with the experimental value? Author's Reply The authors would l ike to thank Prof. Stern for his valuable discussion. F. The computed time period of cavitation cloud shedding, which is nondimensionalized based on the uniform f low velocity and the chord length of the hydrofoi l, is about 1.4 at a c a v i t a t i o n n u m b e r o f 1 . 0 . T h e 683 dimensionalized shedding period is 1.4x(chord length)/(uniform flow velocity). It becomes 0.0117sec (85.7Hz) when the chord length and the uniform flow velocity are 50 mm and 6.0m/s, respectively (the experimental condition). The authors did not measure the time period of the cavitation cloud in their experiments. However, the authors' similar experiment [ Al ] showed that the nondimensionalized time period of cloud cavitation shedding was 2.1. This value is very close to the present computed result. The main purpose of this study was to clarify the generation mechanism of cloud cavitation. The computed result showed a close relationship between the behavior of the separated shear layer and the cavitation cloud. This is because the cavitation cloud shedding frequency also depends on the Reynolds number. Further theoretical investigation is needed to predict the cavitation cloud shedding frequency for arbitrary f low conditions. [Al ] Kubota,A. et al.: Unsteady Structure Measurement of Cloud Cavitation on a Foi l Section Using Conditional Sampling Technique, ASME J. Fluid Eng., Vol.111, No.2, 1989, pp.204-210
