A pressure-based algorithm for cavitating flow computations

42

2011,23(1):42-47 DOI: 10.1016/S1001-6058(10)60086-8

A PRESSURE-BASED COMPUTATIONS*

ALGORITHM

FOR

CAVITATING

FLOW

ZHANG Ling-xin, ZHAO Wei-guo, SHAO Xue-ming Department of Mechanics, State Key Laboratory of Fluid Power Transmission and Control, Zhejiang University, Hangzhou 310027, China, E-mail: [email protected]

(Received September 25, 2010, Revised December 20, 2010) Abstract: A pressure-based algorithm for the prediction of cavitating flows is presented. The algorithm employs a set of equations including the Navier-Stokes equations and a cavitation model explaining the phase change between liquid and vapor. A pressure-based method is used to construct the algorithm and the coupling between pressure and velocity is considered. The pressure correction equation is derived from a new continuity equation which employs a source term related to phase change rate instead of the material derivative of density Dρ / Dt . This pressure-based algorithm allows for the computation of steady or unsteady, 2-D or 3-D cavitating flows. Two 2-D cases, flows around a flat-nose cylinder and around a NACA0015 hydrofoil, are simulated respectively, and the periodic cavitation behaviors associated with the re-entrant jets are captured. This algorithm shows good capability of computating time-dependent cavitating flows. Key words: cavitating flow, pressure-based algorithm, cavitation model, numerical simulation

Introduction Cavitation is often observed in a variety of hydrodynamic systems, such as hydrofoils, marine propellers and underwater bodies. In spite of many efforts, the behavior of cavitation is not yet fully understood. In recent years, large scale computation becomes possible and modeling cavitating flows is of great interest. There are two models usually used for modeling cavitating flows: single-phase model and two-phase model. In the single-phase model, the change in density satisfies a state equation of the single-fluid flow. The two-phase model treats the cavitating fluid as a mixture of water and vapor, capturing cavity with void fraction. The two-phase model is better to simulate complicated cavitating flows due to its more physical considerations[1-4]. Several two-phase void * Project supported by the National Natural Science Foundation of China (Grant No. 10602030), the National Key Basic Research Program of China (973 Program, Grant No. 2009CB724303) and the Fundamental Research Funds for the Central Universities (Grant No. 2010QNA4015). Biography: ZHANG Ling-xin (1978-), Male, Ph. D., Associate Professor

fraction models have been proposed, such as Merkle’s volume fraction model, Singhal’s mass fraction model and Kunz’s volume fraction model[1]. These void fraction models have been widely used to investigate steady and unsteady cavitating flows[5-9]. In this article, Kunz’s volume fraction model is employed. Much attention has been paid to numerical simulation of cavitating flows[10-14]. The major numerical challenge arises due to large density change and its induced compressibility of the mixture. Kubota et al.[15] proposed a Poisson equation for pressure in cavitating flows. However, the expression of this pressure equation was complicated due to their bubble two-phase model. An implicit preconditioned algorithm has been given by Kunz et al.[1] for the computation of viscous two-phase cavitating flows. Senocak and Shyy[16] proposed a pressure-based algorithm, in which the coupling between velocity, pressure and density for proper formulation of the pressure equation was discussed. Although efforts have been made, the methodology of simulating cavitating flows still need to develop. The aim of this article is to propose a new numerical methodology of pressure-based algorithm, which includes the consideration of compressibility

43

due to the phase change. Different from the algorithm proposed by Senocak and Shyy, our algorithm treats pressure coupled with the material derivative of density, instead of the direct coupling between pressure and density. To demonstrate the capabilities of our algorithm, two cavitating flows are simulated and discussed. 1. Governing equations The governing equations include the standard Navier-Stokes equations and a two-phase cavitation model, which is added to take into account the phase change between liquid and vapor. Since there is not a distinct interface between liquid and vapor, a void fraction is introduced to illustrate the evolution of cavitation. Kunz’s volume fraction model[1], based on a transfer equation, is written as follows: ∂ ( ρl α l ) ∂t

+ ∇i( ρlα l u ) = m = m + + m −

(1)

or ∂ ( ρ vα v ) ∂t

+ ∇i( ρvα v u ) = − m

(2)

where the source term m represents the phase change rate, and α l and α v are liquid volume fraction and vapor volume fraction, respectively, which satisfy

αl + αv = 1

(3)

The density of mixture is defined as

ρ m = α l ρl + α v ρv

(4)

In Kunz’s model, the creation of vapor occurs when the pressure is below the saturation pressure and m − is proportional to liquid volume fraction. The destruction rate of vapor, m + , is proportional to vapor volume fraction, but at the same time is a second-order polynomial function of liquid volume fraction. They are given by m+ =

Cprod ρvα l2 (1 − α l ) t∞

(5a)

Cdest ρvα l min[0, p − pv ] (5b) 1 ρlU ∞2 t∞ 2 where Cdest and Cprod are empirical constants, U ∞ m− =

and t∞ are the characteristic variables of the bulk flow. To simulate the turbulence, the two-equation model is used. Thus, the partial differential Reynolds-averaged equations governing unsteady cavitating flows in the Cartesian coordinates may be written in the following general form ∂ ( ρ mφ ) ∂t

+ ∇i( ρ m uφ ) = Sφ

(6)

where all symbols “ − ” for averaged variables are neglected for the ease of presentation, φ is a general variable, representing 1, u1 , u2 , u3 , k , ε , f v , respectively. Here f v is vapor mass fraction, and the relationship between mass fraction and volume fraction is

ρ fv = ρvα v

(7)

Although Kunz’s model is based on liquid volume fraction, a conservative form equation based on vapor mass fraction is obtained and solved together with the Navier-Stokes equations. From Eq.(2) we can get S f = −m

(8)

2. Numerical methodology Generally speaking, water and vapor can be both regarded as incompressible fluids when the bulk velocity is relatively low. However the compressibility of the mixture of water and vapor is complicated due to the special relationship between density and pressure in cavitating flows. As local pressure increases or decreases at the vicinity of saturation pressure, the density of the mixture changes dramatically. The compressibility of cavitating flows can be evaluated by using a state equation of the mixture[17]. In view of the coexistence of the strongly compressible two-phase medium and the incompressible water flow, assuming the whole flow field as either incompressible or compressible may not be a good idea. To cope with this challenge, we propose a new pressure-based algorithm for simulating cavitating flows. We begin with getting a new continuity equation for the cavitating flows, in which the mass change rate is related to phase change rate instead of Dρ / Dt . To do this, Eqs.(1) and (2) can be rewritten as

∂α l m + ∇i(α l u ) = ∂t ρl

(9)

44

ui = ui∗ + ui′

∂α v m + ∇i(α v u ) = − ρv ∂t

(10) After getting p ′ , the estimated

The density of each mass is constant in cavitating flows. For α l + α v = 1 , adding Eq.(9) and Eq.(10) yields

(

)

∇ iu = ρl − ρv m = Cρ S f −1

−1

(11)

∗ 1

where Cρ = ρv−1 − ρl−1

(12)

The algorithm is discussed in the Cartesian coordinates for the ease of presentation. The finite volume method is adopted to discretize the governing partial differential equations. All the unknown variables are located at the cell center node P , and the neighbors of the node P are labeled E , W , N , S , T , B , respectively. The discretized form of the governing equations can be written as

∫ΔV

∂t

dV + [ F1 ] + [ F2 ] + [ F3 ] = e w

t b

n s

∫ΔV Sφ dV (13)

where Fi is the flux of φ at each control volume face, and [ i ] f1 = ( i ) f − ( i f

2

1

)f

2

.

The continuity equation is discretized as

[ M1 ]

e w

+ [M 2 ] + [M 3 ] = n s

t b

∫ΔV Cρ S dV f

(14)

nodes and ui′ at faces

(u ′ )

f

= ( Di ∇i p′ ) f

(18)

where i = 1, 2,3 , respectively, and f denotes each face. Substituting Eq.(17) into the continuity equation, Eq.(14), we obtain e

n

⎡⎣ D1 ( ∇1 p ′ ) a1 ⎤⎦ + ⎡⎣ D2 ( ∇ 2 p ′ ) a2 ⎤⎦ + w s t

⎣⎡ D3 ( ∇3 p′ ) a3 ⎦⎤ b =

∫ΔV Cρ S dV − f

⎧ ⎡ ∗ ⎤e ⎡ ∗ ⎤ n ⎡ ∗ ⎤t ⎫ ⎨ ⎣ M1 ⎦ + ⎣ M 2 ⎦ + ⎣ M 3 ⎦ ⎬ w s b⎭ ⎩

(19)

The first term on the right hand side of the above equation also includes the function of pressure, which plays an important role in the process of phase change. We shall deal with this term in the pressure correction equation. The term Cρ S f can be written as

)

Cρ S f = ξ ⎡ p∗ − pv + p′⎤ + η ⎣ ⎦

(20)

where ξ and η are

ξ =−

(15)

where a1 , a2 and a3 are areas of the cell faces perpendicular to the directions of x , y , z , respectively. Following the procedure employed in the SIMPLE algorithm, the estimated velocity field (u1∗ , u2∗ , u3∗ ) can be calculated out with an estimated

and

∗ 3

(

In the Cartesian coordinates, the expression for M i is

M i = ui ai

∗ 2

p∗

(u , u , u ) can be corrected using the above equations. According to the momentum equations, we have the approximate relationship between p ′ at

i

∂ ( ρφ )

(17)

η=−

Cρ Cdest ρvα l if p < pv 1 2 ρ lU ∞ t∞ 2 Cρ Cprod ρvα l2 (1 − α l ) t∞

(21a)

(21b)

pressure field p∗ . Generally speaking, the estimated velocity components can not satisfy the continuity equation, so the corrections for pressure and velocity are performed as follows:

Note that the expressions of ξ and η depend on the cavitation model, so here in Eq.(21) Kunz’s model is used. Substituting Eq.(20) into Eq.(19) leads to a pressure correction equation. More importantly, the term ξ i p in Eq.(21) should be treated implicitly to improve the stability of the computation. The final pressure correction equation is given by

p = p∗ + p′

APp pP′ = AEp pE′ + AWp pW ′ + ANp pN ′ + ASp pS ′ +

(16)

45

ATp pT ′ + ABp pB′ + SUp

(22)

where APp = AEp + AWp + ANp + ASp + ATp + ABp − ξ iΔV

(23)

e n t ⎫ ⎧ SUp = − ⎨ ⎡ M 1∗ ⎤ + ⎡ M 2∗ ⎤ + ⎡ M 3∗ ⎤ ⎬ + ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ w s b⎭ ⎩

(

)

⎡ξ p∗ − p + η ⎤ iΔV v ⎦ ⎣

(24)

Fig.1 Computational domain and c-grid system around axisymmetric flat-nose cylinder

3. Simulation and discussion Two 2-D examples, an axisymmetric case and a planar case, are simulated for the validation of the algorithm. Time-dependent computations are performed to demonstrate the unsteadiness associated with re-entrant jets and periodic shedding phenomena. An important dimensionless number σ , the cavitation number, is typically defined as

σ=

p∞ − pv 1 ρ U ∞2 2 l

(25)

The axisymmetric case is the flow over a flat-nose cylinder. The cavitating flows around cylinders of various noses have been investigated experimentally by Rouse and McNown, who measured the averaged pressure distributions along cylinder bodies[16]. In our computation, the Reynolds number coincides with that in the experiment, i.e., Re = 1.36 × 105 , and the RNG k − ε model is employed in calculating the turbulent viscosity. Figure 1 shows the configuration and the single-block structural grid (302×182). The free-stream velocity is 10 m/s, and the cavitation number equal to 0.5. A

Fig.2 Transient density contours around flat-nose cylinder at

σ = 0.5

time-sequence of predicted density contours for a flat-nose cylinder simulation is shown in Fig.2, where the white color regions denote the cavities. It can be observed that the entire cavity is highly unsteady, illustrating a periodic evolution with time development. A transport of liquid jet towards the front of the cavity is captured. The moving liquid jet is sheared under the interaction with the vapor region until it reaches the leading edge. The calculated Strouhal number based on the free-stream velocity and cavity length is about 0.3. The computation is performed over the range of 36 periods and the time averaged pressure distribution along cylinder body is illustrated in Fig.3. In the head region of the cylinder, there is a high pressure distribution. A low pressure region appears from the leading edge, which is associated with the development of the cavitation. In the region of cavity closure, the pressure increases moderately. The computed pressure coefficients are well compared with the measured data.

46

in the Strouhal numbers is understandable.

Fig.3 Comparison of pressure coefficient distribution for the flat-nose cylinder at σ = 0.5 , where experimental data are from Ref.[16]

Fig.4 Computational domain and grid for NACA0015 hydrofoil at α = 8

The second case is the cavitating flows over a 2-D NACA0015 hydrofoil. This geometry has been investigated by Kubota et al.[15] and Saito et al. [18]. In our computation, the angle of attack is fixed at α = 8o and the Reynolds number, based on the uniform flow velocity and chord length of the hydrofoil, is 1.0×106. Figure 4 shows the computational domain, the shape of hydrofoil and the computation grid (260×150). The computational unsteady motion of attached cavity of NACA0015 at σ = 1.2 is illustrated in Fig.5. The periodic behaviors also appear, with a clear liquid jet traveling from the end of the cavity up to the leading edge. While the attached cavity shedding, a vortex cavitation, called cloud cavitation, develops downstream, and a new attached cavity grows from leading edge. This periodic behavior is widely observed in experiments. The calculated lift coefficient, shown in Fig.6, illustrated the periodic feature of the hydrofoil cavitation. The present averaged Ct is about 0.56, which is close to the predicted result of 0.54 by Saito et al.[18]. According to the oscillation frequency, the Strouhal number is approximately 0.325. This value is a little higher than Saito’s predicted value of 0.28 for 2-D simulation and their value of 0.3 for 3-D simulation. As it is well known that the cavity length is difficult to determine exactly, the small difference

Fig.5 A periodic evolution of density ( ρ / ρt = 0.9 ) contours around NACA0015 foil at σ = 1.2

4. Conclusions A pressure-based algorithm is presented for turbulent cavitating flow computations. Assume that there is no slip between the two phases, the single-fluid Navier-Stokes equations are employed. That is, the mixture of liquid and vapor is treated as one compressible continuum, and the cavity is captured with void fraction or the density of the mixture. A cavitation model given by Kunz et al. is employed to account for the phase change between two phases.

47

[4] [5]

[6]

[7] Fig.6 Lift coefficient history for NACA0015 foil, with Tc being the characteristic time based on uniform velocity and chord length

General pressure-velocity coupling method is not enough for computing cavitating flows, in which the density changes dramatically. To overcome the difficulty associated with large density change, the relationship between pressure and the material derivative of density is introduced into the continuity equation. The resulting pressure correction equation shares the information of phase change model. Two cavitating flows are computated for the validation of our algorithm. The flow over a flat-nose cylinder is performed at Re = 1.36 × 105 and σ = 0.5 . The computational averaged pressure distribution along cylinder body is well compared with the measured data from Rouse and McNown’s experiment. The re-entrant jet occurring at the end of the cavity is captured. The computation for a NACA0015 hydrofoil flow at the angle of attack of 8o and σ = 1.2 also exhibits the periodic behavior associated with the re-entrant jet. The obtained Strouhal number, 0.325, is close to Saito’s predicted value.

References [1]

[2]

[3]

KUNZ R. F., BOGER D. A. and STINEBRING D. R. et al. A preconditioned Navier-Stokes method for two-phase flows with application to cavitation prediction[J]. Computers and Fluids, 2000, 29(8): 849-875. SINGHAL A. K., ATHAVALE M. M. and LI H. Y. et al. Mathematical basis and validation of the full cavitation model[J]. Journal of Fluids Engineering, 2002, 124(3): 617-624. CHAHINE Georges. L. Numerical simulation of bubble flow interactions[J]. Journal of Hydrodynamics, 2009, 21(3): 316-332.

[8] [9]

[10]

[11] [12]

[13]

[14]

[15]

[16] [17]

[18]

TAMURA Y. MATSUMOTO Y. Improvement of bubble model for cavitating flow simulations[J]. Journal of Hydrodynamics, 2009, 21(1): 41-46. VENKATESWARAN S., LINDAU J. W. and KUNZ R. F. et al. Computation of multiphase mixture flows with compressibility effects[J]. Journal of Computational Physics, 2002, 180(1): 54-77. WU J. Y., UTTURKAR Y. and SHYY W. Assessment of modeling strategies for cavitating flow around a hydrofoil[C]. The Fifth International Symposium on Cavitation. Osaka, Japan, 2003, Cav03-OS-1-007. RHEE S. H., KAWAMURA T. and LI H. Y. Propeller cavitation study using an unstructured grid based Navier-Stokes solver[J]. Journal of Fluid Engineering, 2005, 127(5): 986-994. LINDAU J. W., BOGER D. A. and MEDVITZ R. B. et al. Propeller cavitation breakdown analysis[J]. Journal of Fluids Engineering, 2005, 127(5): 995-1002. PASSANDIDEH-FARD M., ROOHI E. Transient simulations of cavitating flows using a modified Volume-Of-Fluid (VOF) technique[J]. International Journal of Computational Fluid Dynamics, 2008, 22(1): 97-114. LIU De-min, LIU Shu-hong and WU Yu-lin et al. LES numerical simulation of cavitation bubble shedding on ALE 25 and ALE 15 hydrofoils[J]. Journal of Hydrodynamics, 2009, 21(6): 807-813. LI Jie, LU Chuan-jing and HUANG Xuan. Calculation of added mass of a vehicle running with cavity[J]. Journal of Hydrodynamics, 2010, 22(3): 312-318. GONCALVES E., CHAMPAGNAC M. and PATELLA R. F. Comparison of numerical solvers for cavitating flows[J]. International Journal of Computational Fluid Dynamics, 2010, 24(6): 201-216. LI Xiang-bin WANG Guo-yu and YU Zhi-yi et al. Multiphase fluid dynamics and transport processes of low capillary number cavitating flows[J]. Acta Mechanica Sinica, 2009, 25(2): 161-172. BARRE S., ROLLAND J. and BOITEL G. et al. Experiments and modeling of cavitating flows in venture: Attached sheet cavitation[J]. European Journal of Mechanics B/Fluids, 2009, 28(3): 444-464. KUBOTA A., KATO H. and YAMAGUCHI H. A new modeling of cavitating flows: A numerical study of unsteady cavitation on a hydrofoil section[J]. Journal of Fluid Mechanics, 1992, 240: 59-96. SENOCAK I., SHYY W. A pressure-based method for turbulent cavitating flow computations[J]. Journal of Computational Physics, 2002, 176(2): 363-383. IGA Y., NOHMI M. and GOTO A. et al. Numerical study of sheet cavitation breakoff phenomenon on a cascade hydrofoil[J]. Journal of Fluids Engineering, 2003, 125(4): 643-651. SAITO Y., TAKAMI R. and NAKAMORI I. et al. Numerical analysis of unsteady behavior of cloud cavitation around a NACA0015 foil[J]. Computational Mechanics, 2007, 40(1): 85-96.