- Email: [email protected]

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

2010, 22(5), supplement: 893-898 DOI: 10.1016/S1001-6058(10)60048-0

Numerical study of the cavitating flows over underwater vehicle with large angle of attack* Ying Chen 1*, Chuan-jing Lu 1, 2, Jian-hong Guo 1 1

Department of Engineering Mechanics, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiaotong University, Shanghai, China 2 State Key Laboratory of Ocean Engineering, Shanghai Jiaotong University, Shanghai, China * E-mail: [email protected]

ABSTRACT: A Homogenous-Equilibrium-Model based cavitation code was developed to numerically analyze the three-dimensional cavitating flows over underwater vehicle navigating with large angle of attack. RayleighPlesset equation based cavitation model was used together with non-linear eddy-viscosity turbulence model. The computed cavity shapes and pressure distributions along the vehicle body were found to generally accord with experimental results at different cavitation numbers and angles of attack. The forces acting on the vehicle body was studied qualitatively to explain why the body may get damaged during navigation. It was also discovered interestingly that, the variation trend of dragforce along with cavitation number at the conditions of large angle was completely opposite to that at zero angle. KEY WORDS: cavitation; numerical; angle of attack; cavity shape; pressure; force coefficient

1 INTRODUCTION High-speed underwater vehicles which deviate from the moving direction may often get damaged by the asymmetrical cavity appearing around the vehicle body. To reveal the mechanism of such damage and research the characteristics of cavitating flows at the condition of large angle of attack, numerical simulation and analysis was carried out. Along with the development of computer capability,

numerical simulation methods for cavitation flows have been improved greatly. Some cavitation models based on the solving of N-S equation have been adopted, most of which belong to the so called Homogenous-Equilibrium-Model (HEM). Merkle et al.[1] put forward a HEM model based on barotropy relation, ρ = ρ ( p ) . This approach was also used by Byeong et al[2] to predict the period phenomenon of internal flows through duct and nozzle, and used by Wu et al[3] to simulate the laminar cavitating flows around 2D hydrofoil, and used by Chen et al[4] to simulate periodic turbulent cavity shedding in Venturi tube or around hydrofoil. Kubota et al[5] proposed a second kind of model which establishes a relation between local void fraction and mixture density, based on the Rayleigh-Plesset equation. Markatos[6] used a third kind of model of solving equation sets of liquid phase and gas phase respectively. The currently most extensively used models are the HEMs based on the transportation equation of phase fraction, called TEM. Kunz et al [7] proposed a such kind of model to simulate partial sheet cavities. Fu et al[8] used this model to research the cavitating flows over revolution bodies. Chen et al[9] used it to simulate partial or super cavitating flows over underwater vehicle. Singhal et al[10] deduced another TEM from the Rayleigh-Plesset bubble dynamics equation. Chen et al[11] adopted this model to the cavitating flows

Project supported by the National Natural Science Foundation of China (Grant No: 10832007) Biography: CHEN Ying (1979-), Male, Ph. D.

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China over simple revolution bodies and validated the model by comparison with experimental results. Based on this model, the authors developed a reliable and robust computer code for the simulations of complicated cavitating flows. The mathematical model and numerical methods adopted will be narrated in the following parts.

∂u ⎞ 2⎛ Rij = − ρ ui′u ′j = − ⎜ ρ k + μt l ⎟ δ ij ∂xl ⎠ 3⎝ ⎛ ∂u ∂u j ⎞ k + μt ⎜ i + ⎟⎟ + 2β 2 f μ μt ( Sik Ωkj − Ωik Skj ) ⎜ ∂x ε ⎝ j ∂xi ⎠ k −2 β3 f μ μt ( Sik' S kj' − 1 3δ ij Slk' Slk' )

The basic approach grounds on the solving of the standard 3D RANS equations, energy equation, turbulence model equations, and the vapor fraction equation in cavitation model. A SIMPLE typed implicit sequential solver based on 3D body-fit multiblock structured mesh was set up. The RANS equations in integral form are listed below:

∂t ∂ − ∂xi

∂ ( ρ u j ui ) ∂x j

⎛ 2 ∂u j ⎜⎜ p + μ 3 ∂x j ⎝

(1)

=

∂ ⎛ ∂ui ∂u j ⎞ μ⎜ + ⎟ ∂x j ⎜⎝ ∂x j ∂xi ⎟⎠

∂ ( ρu j k )

(2)

⎞ ∂Rij ⎟⎟ + ρ gi + ∂x j ⎠

⎞ ⎟⎟ ⎠

⎛ ∂ε ⎜⎜ ( μ + μt σ ε ) ∂t ∂x j ∂ xj ⎝ 2 ⎛ε ⎞ ⎛ε ⎞ +Cε 1 ⎜ ⎟ ⋅ 2 μt Sij Sij − Cε 2 ρ ⎜ ⎟ ⎝k⎠ ⎝ k ⎠

⎞ ⎟⎟ ⎠

∂t ∂x j +2μt Sij Sij − Ck ρε ∂ ( ρε )

+

∂ ( ρu jε )

Singhal et al’s cavitation model based on the equation of vapor mass fraction f v is expressed as: ∂ ( ρ fv )

=

=

∂ ∂x j

∂ ∂x j

+

∂ ( ρu j fv ) = ∂x j 1

⎡2 ⎤2 p −p , 0) ⎥ (1 − f v ) − Ce ρl ρv ⎢ MAX ( v σ ρl ⎣3 ⎦ k

(6)

1

⎛ ∂k ⎜⎜ ( μ + μt σ k ) ∂x j ⎝

+

and rotation tensor respectively.

⎡2 ⎤2 p − pv , 0) ⎥ f v Cc ρl ρl ⎢ MAX ( σ ρl ⎣3 ⎦ k

A quadratic k − ε turbulence model with modification for low Re conditions was introduced:

∂ ( ρk )

Where, Sij and Ωij denote the velocity strain tensor

∂t

∂ρ ∂ + ( ρu j ) = 0 ∂t ∂x j

+

(5)

ε

2 MATHEMATICAL MODEL

∂ ( ρ ui )

851

(3)

(4)

Where, σ k and σ ε indicate Prandtl numbers. Ck , Cε 1 and Cε 2 denote the model constants with modification for low Re conditions. The extra non-linear Reynolds stress in the RANS equation is:

Where, σ denotes the surface tension of fluid; Ce and Cc are empirical parameters. The source terms are deduced from the Rayleigh-Plesset bubble dynamics equation. This cavitation model reveals the nature of cavitation generation and collapse most essentially, which guarantees its convergence and stability when used to small cavitation number conditions. The mixture density is a function of the mass fractions of liquid and vapor (viz. fl and f v ). The density is calculated by solving f v equation coupled with the RANS equations. The ρ − f relationship is:

1

ρ

=

fl

ρl

+

fv

ρv

,

fl + fv = 1

(7)

3 NUMERICAL METHODS

The solver bases on multi-block non-orthogonal structured grid, which is suitable to kinds of complex geometry. Finite Volume Method (FVM) was used to discretize the computation domain. A colocated variable arrangement on the grid was adopted. Full Multigrid Method (FMG) using V-cycles was also included to accelerate the convergence of iteration. A segregated solver was used, and the implicit method was adopted to discretize the equations. Which means, the variable at volume center φP , and its neighbors,

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

852

φnb , are assumed to be unknown in each outer iteration cycle. Under-relaxation factors were used in the outer iteration when sequentially computing the coupled non-linear equations. Because a mixture hypothesis was supposed in the cavitation model to consider the whole fluid field as a continuous region with variable density, a Pressure-Velocity-Density coupling method of SIMPLE type was established to insure the continuity. For the solving of the transient term in any type of Reynolds’ transportation equation, a three time level method was adopted. Namely, after subtracting the continuity equation from the control equation of φ , the transient term can be approximated as follows: ⎛ ∂φ ⎞ ∫V ⎜⎝ ρ ∂t ⎟⎠dV

⇒ ρ

3φ n +1 − 4φ n + φ n −1 ΔV 2Δt

(8)

The convective flux was composed of first-order upwind scheme and the corresponding deferred correction to compromise convergence, stability, and accuracy. The flux combines first-order implicit part and high-order explicit correction part which is calculated using the variables in the last iteration step. TVD High-Order-Convection scheme was also included in the code to capture great gradient if necessary. The flux can be expressed as: K

K ∫ f ( ρ uφ )dS

(9)

The diffusive flux was also mixed of first-order implicit scheme and deferred correction of central scheme. For surface integral, midpoint rule was adopted. That is, the integral was approximated as a product of the integrand at the cell face center and the cell face area. The phase fraction equation in cavitation model can composes simultaneous equations with the continuous equation to eliminate some terms from the original equation. Thus the coefficient matrix of discretized algebraic equation has the advantage of diagonal dominance and good convergence. If Fme denotes the mass flux through the ‘east’ surface, the upwind implicit parts of the convective terms on all the surfaces of one control volume can be summed as:

∑F k =1

k m

n

∑F

k m

k =1

n

n

k =1

k =1

f k = ∑ Fmk ( f UDS )k = ∑ ( Fmk + f P + Fmk − f nb ) (10)

n

n

k =1

k =1

f k = ∑ Fmk ( f P ) k = ∑ ( Fmk + f P + Fmk − f P )

(11)

The positive convective flux terms can be eliminated by subtracting Eq (11) from Eq (10), and this makes the diagonal coefficients of f in the algebraic equation dominant. The Strongly Implicit Procedure (SIP) method was used to solve the discretized algebra equations iteratively. On multi-block grid, in each inner iteration cycle, the equation set was solved block by block. For control volumes besides block interfaces, the values of neighboring blocks were used explicitly. After discretizing the moment equation implicitly, two computation steps were executed in every time step Δt . The first step is to solve the implicitly discretized momentum equation by iteration. Eq (12) denotes the number- m outer iteration, in which, the number- m intermediate velocity field U m is calculated using p m −1 and U m −1 obtained in the last iteration.

implicit

⇒ ⎡( ρ un ΔS ) f φ UDS ⎤ + ⎣ ⎦ explicit γ ⎡⎣( ρ un ΔS ) f φ HOC − ( ρ un ΔS ) f φ UDS ⎤⎦

n

Where Fmk + and Fmk − respectively represents the mass flux flowing out and in the control volume. The continuous equation can be multiplied by f at both sides and then be discretized, thus can be written as:

uim, P =

ui m Qumi −1 − ΣAnb ui , nb

−

APui

1 ⎛ δ p m −1 ⎞ ⎜ ⎟ APui ⎝ δ xi ⎠ P

(12)

Where Qumi −1 denotes the source term minus the ui contribution of pressure field. APui and Anb are the discretized coefficients of the present point and the neighboring points, respectively. V m does not satisfy the continuous equation, therefore, the second step is to solve the pressure-correction equation, and equation of the imbalance Q m is as follows:

(ρ

m −1

− ρ n ) ΔV Δt

k

(

+ ∑ ρ m −1U m ⋅ n f =1

)

f

= Q m

(13)

The discretized continuous equation and correction relation exist, as the following equations:

(ρ

n +1

k

− ρ n ) ΔV Δt + ∑ ( ρ n +1U m ⋅ n ) = 0

ρ n +1 = ρ m −1 + ρ ′ ,

f =1

U m = U m + U ′

f

(14) (15)

Putting Eq (13) to Eq (15) together, omitting the high-

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

853

order quantities, the total correction equation can be rearranged as: ρ P′ ΔV Δt

k ⎡ ⎛ 1 + ∑ ⎢ Cρ U m ⋅ n p′ − ρ m −1 S ⎜ u f =1 ⎢ ⎝ AP ⎣

(

)

⎞ ⎛ δ p′ ⎞ ⎤ ⎟⎜ ⎟ ⎥ = 0 (16) ⎠ ⎝ δ n ⎠ ⎥⎦

Utilizing the upper equation and the compressibility relation, ρ ′ = Cρ p′ , thus the pressure correction p ′ can be computed out. 4 RESULTS AND DISCUSSION

The present research works were carried out about the cavitating flows around a kind of underwater vehicle navigating with large angle of attack. The model and the computational grid are presented in Fig 1.

Fig. 1 Vehicle model and 3D computational grid

In all the computations below, the working fluid is a liquid at temperature of 293K, liquid and vapor densities are 1000kg/m3 and 0.5542kg/m3 respectively. The saturation pressure is 2350Pa, the surface tension is set to be 0.0717N/m, and the incoming velocity U ∞ is adopted as 10m/s (Re ≈ 105). Working condition of σ =0.01~0.6 and α =0°~20° were studied in detail.

Fig. 2 Computational and experimental results at varieties of angles of attack ( σ =0.3)

Fig 2(a) presents the computed 3D iso-surfaces ( f v =0.1) of the cavities over the vehicle at varieties of angles of attack, and for each angle three views are [12] shown. The experimental photos by Quan et al were also shown in Fig 2(b), and in each case two views are shown. By comparison, it is clearly indicated that the computation results have good agreement with the experimental ones in cavity shape. It obvious to see that, the computational results and the experimental images are similar at the following aspects: (1) For identical angle of attack, the cavity lengths on the three sides in turn are, the length on the side direction, on the suction side, and on the pressure side; (2) Cavity lengths on the suction side change little for different angles; (3) When the angle increases, cavity on the side direction prolongs, while the cavity length on the pressure side decreases obviously.

Fig. 3 Pressure distribution for varieties of angles of attack

854

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

Fig 3 presents the simulated pressure distribution along the pressure-side wall and suction-side wall at different cavitation numbers and angles of attack, comparing with the experimental results obtained by Quan et al[12]. The computational results generally agree with the experimental data well on cavity length, pressure inside cavity and pressure at cavity closure point et al. What’s needed to be concerned is, the experimentally measured pressure in all-wet area after the cavity closure, especially on the suction wall at condition of σ =0.3, is confused to some extent. The experimental and the numerical results reveal the following identical pressure distribution characteristics: along the axis of the vehicle, firstly, the pressure reduces intensely from the maximal value at the stagnant point to the saturated vapor pressure p∞ at the forepart separating point of cavity. Then, pressure keeps constant of p∞ throughout the pressure platform area inside cavity. Finally, the pressure lifts dramatically across the cavity closure point and then drop lightly to keep constant value at the following all-wet part of the vehicle wall. The relation between the length of pressure platform section and the angle of attack presents that, the cavity length on the pressure side decreases along with the increasing of the angle of attack. That is said, the cavity length on pressure side decreases along with the angle more apparently than on suction side. This obvious character is also successfully simulated to be in accordance with experiment.

from its moving direction seriously and large angle occurs, the vehicle body may usually bear great load or even get damaged. To explain the mechanism of such phenomenon, a circumference integral stress coefficient CFp was defined as: 2π

CFp

∫ ( p − p )r sin θ dθ = ∞

0

1 2

(17)

ρU ∞2 ⋅ 2πr

Fig 4 presents the distribution of CFp along vehicle axis for conditions of different angles of attack. Contrasting with cavity image and pressure distribution, it’s obvious that the stress peaks appearing in this figure are caused by asymmetry of cavity on pressure side and suction side. As the pressure behind cavity closure point does not differ greatly for one specific cavitation number when the angle decreases, the peak stresses caused by difference between this pressure and saturated vapor pressure are basically similar. However, when the angle rises, the extent of positive stress widens along with strengthening of cavity asymmetry. Table 1 Resultant force on unit length at different conditions Cavitation number 0.2

10 m/s 8×104 N

Incoming velocity 20 m/s 30 m/s 3×105 N 7×105 N

50 m/s 2×106 N

0.3

1×105 N

5×105 N

3×106 N

1×106 N

Taking vehicle radius of r =0.25m as example, Table 1 lists the resultant force on unit length calculated using peak stress. It’s obvious to see from the resultant force that the bending moment acting on vehicle body is so large to influence the vehicle or even damage it.

Fig. 5 Relation between force coefficients and angle of attack at different cavitation numbers

Fig. 4 Resultant stress coefficient along vehicle axis ( σ =0.3)

In real applications, when underwater vehicle deviates

Fig 5 presents the relation of lift-force and drag-force coefficients with the angle of attack at different cavitation numbers, comparing with the relation of allwet flow condition. An important and interesting conclusion was found that, the relation between the drag-force coefficient and cavitation number at the condition of large angle of attack is absolutely

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China negative to the relation at zero angle condition. Such result can be explained as: At zero angle condition, the cavity is axisymmetric, thus drag force is primarily contributed by viscous force. Therefore, the smaller the cavitation number is, the longer the vehicle is surrounded by cavity, thus the smaller the drag force is. However, at large angle conditions, the cavity is asymmetric, thus drag force primarily results from the difference between the pressures on suction side and pressure side. Hence, along with the increment of cavity length on suction side and decrement on pressure side, the pressure difference increases and thus lift force and drag force increase. 5 CONCLUSIONS

The HEM based mathematical models and numerical simulation methods for cavitation were established, and a computation code suitable for the computation of varieties of complicated cavitating flows was developed. The authors applied the methods and computer code to research the 3D cavitating flows over a kind of underwater vehicle navigating with large angle of attack. Full field density, pressure and velocity distributions at different cavitation numbers and different angles of attack were studied systematically. The computed cavity shapes and pressure distributions agree with experiment in general and inflect some detailed cavity characters observed by experiment. The integral stress acting along vehicle axis was studied and found that the resultant stress in real flows may be great enough to damage the vehicle. The mechanism was considered to result from asymmetricy of cavity and explained by pressure difference of suction side and pressure side. It was also discovered interestingly that, the dragforce variation trend along with cavitation number at large angle conditions was absolutely opposite to that at zero angle condition.

855

REFERENCES [1] MERKLE C L, FENG J, BUELOW P E O. Computational modeling of the dynamics of sheet cavitation [C]. 3rd Int. Symp. on Cavitation. Grenoble, France, 1998: 307-313. [2] BYEONG R S, SATORU Y, XIN Y. Application of Preconditioning method to gas-liquid two phase flow computations [J]. J Fluids Engineering, 2004, 126(4): 605612. [3] WU L, LU C J, LI J et al. Numerical simulations of 2D periodic unsteady cavitating flows [J]. Journal of Hydrodynamics, Ser B, 2006, 18(3): 341-344. [4] CHEN Y, LU C J, WU L. Modeling and computation of unsteady turbulent cavitation flows [J]. Journal of Hydrodynamics, Ser B, 2006, 18(5): 559-566. [5] KUBOTA A, KUTO H, YAMAGUCHI H. A new modeling of cavitating flows: a numerical study of unsteady cavitation on a hydrofoil section [J]. J Fluid Mech, 1992, 240(7): 59-96. [6] MARKATOS N. Modeling of two-phase transient flow and combustion of granular propellants [J]. Int J Multiphase Flow, 1986, 12: 913-933. [7] KUNZ R, BOGER D, CHYCZEWSKI T et al. Multiphase CFD analysis of natural and ventilated cavitation about submerged bodies [C]. 3rd ASME/JSME Joint Fluids Engineering Conference, San Francisco, CA, USA, 1999. [8] FU H P, LU C J, WU L. Research on characteristics of flow around cavitating body of revolution [J]. Journal of Hydrodynamics, Ser A, 2005, 20(1): 84-89. (in Chinese) [9] CHEN Y, LU C J, Wu L. Numerical method for three dimensional cavitating flow of small cavitation number [J]. Chinese Journal of Computational Physics, 2008, 25(2): 163-171. (in Chinese) [10] Singhal A K, Athavale M M, Li H Y et al. Mathematical basis and validation of the full cavitation model [J]. J Fluid Engineering, 2002, 124: 617-624. [11] CHEN Y, LU C J. A Homogenous-Equilibrium-Model based numerical code for cavitation flows and evaluation by computation cases [J]. Journal of Hydrodynamics, Ser B, 2008, 20(2): 186-194. [12] QUAN X B, LI Y, WEI H P et al. An experiment study on cavitation of underwater vehicle’s surface at large angles of attack [J]. Journal of Hydrodynamics, Ser A, 2008, 23(6): 662-667. (in Chinese)