Towards numerical prediction of unsteady sheet cavitation on hydrofoils

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China

699

2010, 22(5), supplement :741-746 DOI: 10.1016/S1001-6058(10)60024-8

Towards numerical prediction of unsteady sheet cavitation on hydrofoils Da-qing Li *, Mikael Grekula, Per Lindell SSPA Sweden AB, Göteborg , Sweden * E-mail: [email protected] ABSTRACT: The paper presents a study of using a modified k-ω model to predict the unsteady cavitating flows around 2D and 3D hydrofoils in the framework of multi-phase mixture flow RANS approach. The cavitation is modeled by SchnerrSauer’s cavitation model. A 2D NACA0015 foil at cavitation number σ=1.0 (unsteady with cloud shedding) is studied first, followed by the Delft twisted 3D foil. It is found that the present RANS method is able to predict the essential features like re-entrant jets, the periodic shearing and shedding of cloud cavities. Two distinct shedding dynamics are noted for the 2D foil: (a) Shedding of medium to large scale structures (at low frequency); (b) Shedding of a series of secondary vortex cavities (at high frequency). For the 3D twisted foil, the collaborated effect of re-entrant jets from the curved closure line to break up a primary cavity, as well as the formation, rollup and transport of cavitation vortices that are observed in the experiment are truly reproduced in the simulation. The method is found to have a tendency to under-predict the lift coefficients. KEY WORDS: cloud cavitation; shedding; turbulence model; cavitation model; hydrofoil.

1 INTRODUCTION Fluctuating sheet cavitation on marine propellers and hydrofoils calls for special concern as it not only degrades hydrodynamic performance but also induces flow instability, vibration and noise. Furthermore, if the sheet cavity itself or its shed structures are performing a very fast collapse at a focal location near a foil surface, cavitation erosion may take place on the hydrofoil [1]. Although experimental observation and measurement using advanced techniques are still indispensable for studying and solving cavitation problems today, advancement of computer power and numerical methods open up new possibilities to study the behavior of cavitation dynamics in greater detail than what is ever affordable by experiments. A brief review of the experimental and computational work for the sheet cavitation was given in reference [11]. The present work is a numerical study using a singlefluid RANS method based on the assumption of homogenous equilibrium mixture of multi-phases

flows. Several recent work and our own studies reveal that (a) the cavitating flow in the mixture-phase region is locally compressible [7][9][10]; (b) the standard twoequation turbulence models (e.g. k-ε class) have a tendency to over-predict the turbulent viscosity at the closure of cavity and its downstream region in reference[3],11]. It has been suspected that the inherent unsteadiness of cavitation is significantly dampened by the viscosity and the models are unable to capture the shedding dynamics. The deficiency of standard models were also observed for example in [5][9][15]. In our previous study [11], a modification was implemented in a RNG k-ε model following the idea of Reboud et al. [5]. This led to a fairly good prediction of re-entrant jets and shedding of cloud cavities on the 3D twisted Delft hydrofoil when compared with the experiment results [2]. The study manifests that, in addition to the physical modeling of cavitation process itself, the turbulence modeling has a critical role in the cavitation prediction. For non-cavitating turbulent flows, the k-ω turbulence models are known to have better performance over k-ε models for wall-bounded boundary layer flows under adverse pressure gradient and free shear flows [14]. It is of engineering interest and significance to study their capability to predict cavitating flows and the performance breakdown due to cavitation. The aim of the work is thus to study the feasibility to use a modified k-ω model (of Wilcox) associated with Schnerr-Sauer’s cavitation model [13] to predict the unsteady cavitation behavior of hydrofoils. The mixture flow RANS solver in FLUENT12 is used in the study. 2 NUMERICAL MODELS 2.1 Multi-phase mixture model The mixture model assumes the working medium

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China

700

being one-fluid with homogeneous mixture of two phases (liquid and vapor). Therefore only one set of RANS equations is needed for the mixture fluid. Denoting the density of the mixture fluid by ρm, the continuity equation for the mixture flow becomes: ∂ G (1) ( ρ m ) + ∇ ⋅ ( ρ m vm ) = 0 ∂t The momentum equation for the mixture reads: G G G ∂ ( ρ m vm ) + ∇ ⋅ ( ρ m vm vm ) ∂t (2) G GT G ⎡ ⎤ = −∇p + ∇ ⋅ ⎣ μ m ( ∇vm + ∇v m ) ⎦ + ρ m g The mixture density is composed of: ρ m =α v ρ v + (1 − α v ) ρl (3) where αv and αl are the volume fraction of vapor and liquid respectively. To close the equation system an additional transport equation is solved for αv. ∂ G (α v ρ v ) + ∇ ⋅ (α v ρ v vm ) = Re − Rc (4) ∂t where the source terms Re and Rc are used to account for the mass transfer between phases and need to be modeled via a cavitation model. 2.2 Cavitation model

The cavitation model proposed by Schnerr and Sauer [13] defines the source term Re and Rc as follows, ρρ 3 2 ( pv − p ) Re = v l α (1 − α ) (5) ρm ℜB 3 ρl when pv > p, and ρρ 3 2 ( p − pv ) Rc = v l α (1 − α ) (6) ρm ℜB 3 ρl when pv < p. The bubble radius RB can be determined by: 1

⎛ α 3 1⎞3 ℜB = ⎜ (7) ⎟ ⎝ 1 − α 4π n0 ⎠ where n0 is the bubble number density given as a constant. The default value n0 = 1013 is used.

used in Eq. 9. The QUICK scheme is used for convection terms in RANS equations. The time step is 5e-4 s and 20 iterations are run within each time step. 3 GRIDS AND FLOW CONDITIONS 3.1 2D NACA 0015 foil The foil has a 0.2m chordlength(C) and is located in the middle of a 0.57m high rectangular tunnel. The inlet plane of the computational domain is situated at 2.0C upstream the leading edge (L.E.) and the outlet plane at 4.0C downstream trailing edge (T.E.). It is operated at 6° angle of attack (AoA) and a cavitation number σ=1.0 Table 1 gives the flow and boundary conditions. The physical properties is from a water temperature at 24°C. No-slip boundary condition is specified on foil surfaces whereas slip condition is assumed on tunnel walls for all cases. The grid sensitivity was studied by three similar grids in our previous work [3]. It was found in general that the finer the mesh, the greater detail of the cavity and vortex structures could be captured. The vapor-liquid interface becomes a bit sharper with finer grids. But the difference in grid resolution does not change the underlying flow pattern and shedding behavior. The total lift and drag predicted on three grids in the noncavitating conditions are almost the same [3]. Due to time constraint, only the results on the coarse grid are available at the time of writing. Fig. 1 shows the coarse grid consisting of 35030 cells. It is evident here that even the coarse grid is substantially refined. Table 1. Flow and boundary conditions for 2D NACA 0015 Boundary Conditions Values Inlet: constant velocity, [m/s] 6.0 Outlet: pressure outlet, [kPa] 20.9 Physical properties of water Vapour Liquid Vaporization pressure, [kPa] 2.97 Dynamic viscosity, [kg/ms] 9.95x10-6 1.10x10-3 3 Density, [kg/m ] 0.023 998.0

2.3 Turbulence model

The formulation and model constants for Wilcox’s k and ω equations are unchanged in the present method. Instead, a modification is made to the formula for turbulent viscosity μt following the idea of CoutierDelgosha and Reboud et al. [7]. k μt = f ( ρ )α * , (8)

ω ( ρ − ρ v )n , where n > 1 f (ρ ) = ρv + m ( ρ l − ρ v )n −1

(9)

The α* is the original damping coefficient. The function f(ρ) aims to rectify the otherwise overpredicted turbulent viscosity in the mixture region by the standard model. An exponent constant of n=100 is

Fig. 1 Coarse grid for NACA0015 foil

3.2 3D delft twisted hydrofoil

The twisted foil is a rectangular wing (0.15m chord by 0.3m span) of NACA0009 profile, Fig. 2(a). The foil is twisted spanwise to give a maximum AoA= 9° at mid-span section and AoA= -2° at the tunnel wall. The cavitation will develop mainly on the suction side mid-span area. The foil is operated at σ=1.07 with a freestream velocity 6.97m/s. The tunnel has a 0.3m by 0.3m cross section. The inlet plane is located at 2.0C

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China upstream the L.E. and the outlet plane at 4.0C downstream the T.E., Fig. 2(b). Table 2 gives the flow and boundary conditions. An O-H type grid is generated for the domain with sufficient refinement (y+≈1) towards foil surface. The grid has about 1.04 million cells.

701

field for a sequence corresponding to Fig. 3 (d).

Fig. 2 (a) Geometry of 3D twisted hydrofoil; (b) Surface mesh of the computational domain Table 2. Flow and boundary conditions for 3D twist foil Boundary Conditions Values Inlet: constant velocity inlet, [m/s] 6.97 Outlet: pressure outlet, [kPa] 29.0

4 RESULTS AND DISCUSSIONS 4.1 2D NACA 0015 foil at σ=1.0

Despite that the foil is operated in a steady inflow, the cavitating flow around the foil at this condition is highly unstable, exhibiting two different shedding phenomena. (1) Intermittent shedding of medium to large scale cloud cavities; (2) Shedding of a series of small cavitating vortices downstream an upstream moving sheet cavity, to be discussed below. 4.1.1 Shedding of medium to large structures

Fig. 3 shows 8 sequences of a shedding cycle. The process starts with an attached sheet developed extensively up to 65%C. It then sheds a first cavity at its closure without the presence of re-entrant jet until Fig. 3 (b). The re-entrant jet is seen to form below the tail of sheet cavity in Fig. 3(c) and starts to move upstream. The shearing effect created by the re-entrant jet and the external main flow breaks off a large piece of cavity at 40%C, Fig. 3 (d). The re-entrant jet appears to be thin, strong and quite durable so that it continues going upstream and breaks off the front part of remaining sheet, Fig. 3 (e)-(f). Meanwhile, the shed cavity rolls up into a large volume and becomes cloudy as it is transported further downstream and leaves the foil, Fig. 3 (g)-(h). This kind of the shedding is mainly caused by the upstream moving reentrant jet and its interference with the external flow, i.e. an effect called “shear by filling the sheet with internal jets” according to Bark et al.[1]. To show the presence of re-entrant jet, Fig. 4 depicts the velocity

Fig. 3 Contours of vapor void fraction showing cavity shape and shedding sequence. Red color= vapor, blue= liquid

Fig. 4 Velocity field of sequence (d) in Fig. 3. Color code: blue for pure liquid; red for pure vapor

Four photos from an experiment [17], performed at AoA=8° and σ=1.2, are shown in Fig. 5 to highlight qualitatively the similarity of experimentally observed shedding pattern with the predicted one.

Fig. 5 Experiment on NACA 0015 at AoA=8° and σ=1.2 [17]

4.1.2 Shedding of a series of cavitating vortices

The process starts with development of a partial cavity with an “open” closure region and without any noticeable re-entrant jet. Due to possibly downstream disturbance of previously shed cavities, the attached sheet performs a fast upstream collapse (towards leading edge) in the same time it releases a series of small vortex cavities, the sequence of which is shown in Fig. 6. The last experiment photo in Fig. 5 seems to

702

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China

support this shedding phenomenon as well.

Fig. 6 Shedding of cavitating vortices for an upstream moving collapse sheet. Color code: red =vapor, blue = liquid

The time history of vapor volume is shown in Fig. 7. It is interesting to note that the large volumes at long periods (low frequency) correspond closely to the shedding phenomenon (1) whereas the small vapor volumes of short periods (high frequency) correlate with the phenomenon (2) described above. The FFT analysis of lift oscillations yields a multiple shedding frequency shown in Fig. 8. Noticeable frequencies include 2.5, 4.6, 8.0 and 14.5Hz, with dominance at low frequency region (<5Hz). This is contrary to the commonly cited one shedding frequency (≈16Hz) in literature. Mesh refinement in another study [11] did not change the situation. However, this is not the case for 3D twisted foil where only one frequency is dominant (to be shown below). The reason for such a difference in frequency for the 2D foil is yet to be clarified. The time-averaged lift (CL) and drag coefficient (CD) are 0.499 and 0.0837, showing a difference by -16% for CL and 14% for CD if compared with the averaged data over AoA=5° and 7° in [16]. For the noncavitating case on the same mesh, the accuracy of the computed CL is well within 2.5% of experiment data [11] . Although one may argue that the high oscillation of lift in cavitating conditions makes the uncertainty of measurement higher, yet the relatively large difference from the measurement indicates that the lift coefficient CL is under-predicted.

Fig. 7 Variation of total vapor volume with time, 2D foil.

Fig. 8 Frequency spectrum based on the lift coefficient CL

4.2 3D delft twisted foil Due to the twist design and the maximum AoA at mid-span, sheet cavitation develops mainly in the mid-span area near the leading edge and has a curved closure line. As a result, the re-entrant jet is no longer purely a reversed flow going upstream but has also a transverse component in spanwise direction. The collaborated effect of these two components causes a highly complex shedding. 4.2.1 Shedding dynamics Fig. 9 compares the shedding dynamics of the predicted cavity (right) with the images (left) of the high-speed video [2] in one cycle. The predicted cavity interface is represented by an iso-surface of vapor volume fraction αv=0.1. The streamlines are plotted to visualize various flow directions. The shedding can be broken up into three events with some minor overlaps: (1) Upstream interface disturbance, Fig. 9 (a)-(b) Fig. 9 (a) shows a state where the sheet cavity has reached approximately its maximum length. The upstream-going re-entrant jet approaches the cavity detachment line and has created many small bubbly disturbances on the cavity surface, yet the sheet remains attached as a whole. The small bubbly disturbances are not captured by the CFD method. Instead, it shows that some large parts are already pinched off by the re-entrant jet. The streamlines clearly visualize the existence of re-entrant jet and the jet front location. What can be also seen from (a) and (b) is a small isolated cloud cavity near the trailing edge, shed from the previous cycle. (2) Primary shedding, Fig. 9 (c)-(f) In Fig. 9 (c)-(e), a large (primary) cavity is cut off and separated from the remaining attached sheet. The CFD shows a very similar pattern. (3) Secondary shedding, roll-up, convection, (g)-(j) In Fig. 9(g), the remaining sheet starts to growth again. The primary shed cavity is convected downstream and starts to roll up and becomes cloudy due to change of pressure. The shed cavity predicted by CFD is slightly smaller but has the same behavior, Fig. 9(h). At the same time a secondary shedding occurs at the closure of the attached cavity – a pair of small cavities are shed, Fig. 9(i). For CFD in Fig. 9(j), the missing parts at the tail of sheet cavity imply that the secondary shedding has taken place, but the shed cavities themselves are not seen. Another difference is that the primary shed cavity starts to collapse earlier in the simulation than in the experiment, therefore the cloud cavity looks smaller. Fig. 10 compares the PIV images in experiment [2] with the predicted cavity at the mid-span section plane in one cycle. As seen from the sequences, the predicted cavities are in general somewhat smaller.

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China

703

Fig. 10 PIV image (left) vs. vapor void fraction at mid-span section plane, 3D twisted foil

Particularly, the predicted cloud cavity does not survive as far downstream as shown in the PIV image, indicating a likely difference in the pressure field between the computation and the test. The development of cavity, shedding and roll-up are very dynamic and none of the cycles are identical, which is true for the experiment observations as well as the simulations. Nevertheless, the periodic behavior is well predicted using the present method. 4.2.2 Global variables

Fig. 9 Top view of one shedding cycle. Left: video image. Right: iso-surface αv=0.1. Flow is from top to bottom.

The time-averaged CL is 0.43, about 15% lower than the measured value of 0.51 [2]. The RANS method shows a tendency to under-predict the lift also for this 3D foil. The pressure was measured at limited number of points on suction side. Fig. 11 gives an example of the time-averaged Cp distribution on the mid-span section compared with the measured data. As pointed out in the preceding paragraph, the predicted pressure downstream the sheet cavity may be different from the experiment one. The difference is seen here in the region between 0.4C and 0.6C where the computed pressure (red line) recovers earlier and quicker than the data, thus leading to an earlier collapse of the primary shed cavity. The time history and frequency spectrum of the total vapor volume are shown in Fig. 12(a) and (b). The dominant shedding frequency is found to be 25.1Hz, comparable to the experimental value of ~21Hz.

704

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China in predicting performance breakdown of hydrofoils. This remains to be a challenge for RANS methods. REFERENCES [1]

[2] Fig. 11 Time averaged Cp distribution at mid-span section [3]

[4] [5]

[6]

Fig. 12 (a) Time history of the total vapor volume, and (b) Spectral analysis of the vapor volume variations

5 CONCLUSIONS

Using a modified k-ω turbulence model combined with Schnerr-Sauer’s cavitation model in a multiphase mixture-flow RANS solver, the shedding dynamics of unsteady partial cavity on a 2D NACA0015 foil and a 3D twisted foil is numerically studied. The results justify that the present RANS method is capable to capture the crucial features like re-entrant jets, periodic shedding and collapse of cavities, as well as the formation, roll-up and transport of cavitating vortices. They are qualitatively in well agreement with experiment observations. As compared with use of a standard turbulence model [3], it has been a significant advancement towards numerical prediction of unstable sheet cavitation. For 2D NACA0015 foil at σ=1.0, multiple shedding frequencies are found in simulation with the most dominant frequency being less than 5Hz. Only the last frequency (≈14.5Hz) is close to the commonly cited frequency (~16Hz) in the literature. For 3D twisted foil, the predicted shedding frequency is 25.1Hz, in reasonable agreement with the experimentally observed value (~21Hz). Common to both foils, the lift coefficients are significantly under-predicted in cavitating conditions (but not in non-cavitating conditions). Further improvements are needed in the cavitation and turbulence model to achieve a higher level of accuracy

[7]

[8]

[9] [10]

[11]

[12] [13] [14] [15]

[16] [17]

Bark G, Grekula M, Bensow R E, et al. On Some Physics to Consider in Numerical Simulation of Erosive Cavitation [C]. The 7th Int. Symp. on Cavitation (CAV2009), Ann Arbor, USA, 2009. Foeth E –J, Terwisga T van. The structure of unsteady cavitation. Part I: Observations of an attached cavity on a 3dimensional hydrofoil [C]. The 6th Int. Symp. on Cavitation (CAV2006), The Netherlands, 2006. Li D -Q, Grekula M, Lindell P. A modified SST kω Turbulence Model to Predict the Steady and Unsteady Sheet Cavitation on 2D and 3D Hydrofoils [C]. The 7th Int. Symp. on Cavitation (CAV2009), Ann Arbor, USA, 2009. Venkateswaran S, Lindau J, Kunz R, et al. Computation of Multiphase Mixture Flows with Compressibility Effects [J]. Journal of Computational Physics, 2002,180: 54-77, Reboud J L, Stutz B, Coutier O. Two-phase flow structure of cavitation: Experiment and modelling of unsteady effects [C]. Proc. of the 3rd Int. Symposium on Cavitation, Grenoble, France, 1998. Coutier-Delgosha O, Fortes-Patella R, Reboud J L. Simulation of unsteady cavitation with a two-equation turbulence model including compressibility effects [J]. Journal of Turbulence, 3 (2002) 058, 2002. Coutier-Delgosha O, Reboud J L, Fortes-Patella R. Evaluation of the Turbulence Model Influence on the Numerical Simulations of Unsteady Cavitation [J]. Journal of Fluids Engineering, 2003, 125: 38-4. Dular M, Bachert R, Stoffel B, et al. Numerical and experimental study of cavitating flow on 2D and 3D hydrofoils [C]. The 5th Int. Symposium on Cavitation (CAV2003), Osaka, Japan, 2003. Wu J, Utturkar Y, Senocak I, et al. Impact of turbulence and compressibility modeling on three-dimensional cavitating flow computations [C]. AIAA 2003-4264, 2003. Wu J, Wang G, Shyy W. Time-dependent turbulent cavitating flow computations with interfacial transport and filter-based models [J]. Int. J. Numer. Meth. Fluids, 2005, 49: 739-761. Li D-Q, Grekula M.. Prediction of dynamic shedding of cloud cavitation on a 3D twisted foil and comparison with experiments [C]. The 27th Symposium of Naval Hydrodynamics, Seoul, Korea, 2008. Huuva T. Large eddy simulation of cavitating and noncavitating flow [M]. PhD thesis, ISBN/ISSN: 9789173850636, Chalmers Univ. of Tech., Sweden, 2008. Schnerr G H, Sauer J. Physical and Numerical Modeling of Unsteady Cavitation Dynamics [C]. In Fourth International Conference on Multiphase Flow, New Orleans, USA, 2001. Wilcox D C. Turbulence Modeling for CFD [M]. DCW Industries, Inc., La Canada, California, 1998. Basuki W, Schnerr G H, Yuan W. Single-phase and Modified Turbulence Models for Simulation of Unsteady Cavitating Flows [C]. Proc. of the German-Japanese Workshop on Multi-Phase Flow, Germany, 2002. Kjeldsen M, Arndt R E A, Effertz M. Spectral Characteristics of Sheet/Cloud Cavitation [J]. J. Fluid Eng., 2000, 122: 481-487. Yakushiji R, Yamaguchi H. Investigation for unsteady cavitation and re-entrant jet on a foil section [J]. J. of the society of naval architects of Japan, 2001, 189.