Numerical modeling of cavitation characteristics and sensitivity curves for reversible hydraulic machinery

Engineering Analysis with Boundary Elements 41 (2014) 18–27

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Numerical modeling of cavitation characteristics and sensitivity curves for reversible hydraulic machinery Anton Iosif, Ioan Sarbu n Department of Building Services Engineering, Polytechnic University Timisoara, Piata Bisericii 4A, 300233 Timisoara, Romania

art ic l e i nf o

a b s t r a c t

Article history: Received 8 April 2013 Accepted 20 December 2013 Available online 28 January 2014

Details on the existence, extent and effects of cavitation can be helpful during the design stages of hydraulic machinery to minimize the effects of cavitation and optimize the design. Fluid flow in reversible hydraulic machinery elements is a complex three-dimensional problem. In this paper, an explicit numerical model based on finite element method (FEM) - dual reciprocity method (DRM) coupling is developed for the 2D simulation of the flow velocity and pressure distributions on the runner blade of a Francis-type reversible radial-axial hydraulic machine. The model is based on the assumptions of an ideal incompressible fluid and relative rotational motion. Instead of directly simulating the 3D flow, the proposed model is based on the analysis of specific 2D flows: fluid axial-symmetric motion and the stream function in an associated plane. This model is applied to reversible radial–axial hydraulic machinery that operates as a pump. The blade exhibits the basic profile NP205. The skeleton of this profile is defined by a quadratic equation, and its thickness function exhibits a NACA profile with a maximum relative thickness of five percent. The numerical results for different discharge values can be used to determine the cavitation characteristics and sensitivity curves for reversible hydraulic machinery. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Francis pump-turbine Radial-axial profile cascades FEM-DRM coupling Flow simulation model Cavitation characteristics and Sensitivity curves

1. Introduction Hydraulic machinery is an essential component of modern technology. Its versatility, high efficiency and minimum environmental impact contribute to its irreplaceable role as a high-quality source of peak power (and energy storage) in electricity generation systems, as well as a versatile power link in hydraulic pressure circuits for diverse technical applications. Reversible machines for extremely high specific energies have been developed due to the demand for total economic and technical optimization; in addition, special machines with increasing capacity have been developed for low-head sites. Reversible hydraulic machinery are commonly applied where low pressures are routinely generated by machine action, e.g., on blade surfaces, with a consequent possibility of cavitation. Cavitation in hydraulic machinery presents unwanted consequences, such as flow instabilities, excessive vibrations, damage to material surfaces and degradation of machine performance. The ability to model cavitating flows has generated significant interest in the CFD community. Details of the existence, extent and effects of cavitation can be helpful during the design stages of hydraulic machinery to minimize the effects of cavitation or to optimize the design. Considerable research on cavitation within pumps and hydropower turbines has been performed over the past decade and extensive reviews are available in the literature [1–5].

n

Corresponding author: Tel.: þ 40 256403991; fax: þ 40 256403987. E-mail address: [email protected] (I. Sarbu).

0955-7997/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.12.007

The fluid motion in reversible hydraulic machinery elements is a complex three-dimensional (3D) problem [6–12]. Computation of 3D flow in rotating and stationary blade passages of turbo-machinery is approximated by assuming that 3D flow can be represented by separate and approximately orthogonal two-dimensional (2D) flows [13,14]. Modern computational techniques facilitate the solution of 2D flows with imposed boundary conditions using different numerical methods [15–18]: the finite element method (FEM), the boundary element method (BEM), and the dual reciprocity method (DRM). According to the authors, coupling the FEM with the DRM enables the transformation of the 3D fluid flow in the reversible hydraulic machinery (pump-turbine) runner to a simpler 2D problem assuming an ideal incompressible fluid. The following stages are recommended: (a) Axisymmetric potential motion, which is solved for a given domain, can be employed to determine the hydrodynamic field as well as the velocity and pressure distributions along the streamlines using FEM. This method prevents the singularities generated by points on the symmetry axis Oz. The FEM algorithm is robust and generally converges with reasonable amounts of computation. (b) The fluid motion around the radial-axial profile cascades, which are disposed on the stream surfaces, is analyzed using the DRM. Unlike the domain meshes required for the FEM, the DRM only discretizes the boundary. By applying this method, onl6y the values of the stream function Ψ and normal derivatives of these functions (tangential velocities) in the nodes, where they are unknown, are obtained. Compared with the FEM, the DRM also

A. Iosif, I. Sarbu / Engineering Analysis with Boundary Elements 41 (2014) 18–27

preserves the advantages of the BEM: a shorter computational time and a reduced number of boundary elements. The following assumptions are considered in solving the 2D flows: – – – – –

fluid is inviscid and incompressible; mass forces are negligible; absolute motion is potential and stationary; rotation number of runners is constant; motion is considered axial-symmetric in the absence of runner blades; – stream surfaces are revolving surfaces; and – motion is uniform at half of the upstream cascade spacing and downstream of the cascade. The relative fluid motion around radial-axial profile cascades is rotational. To determine this motion, the conformal mapping of the domain on a stream surface in an associated plane is employed, as demonstrated by Abdallah [13] and Carte [19]. The crossing of a stream surface with the runner blades generates a radial-axial cascade of profiles. Because a stream surface is deployable into a plane surface, it is conformally mapped onto an associated plane (Prasil), whereas the radial-axial cascade is mapped onto a linear plane. The motion in this plane is determined by solving a boundaryvalue problem for the differential equation of the stream function ψn and obtaining the velocity and pressure field in the associated plane. The results are subsequently transposed to the stream surface. The numerical simulation model is applied to the pump radialaxial cascade that is presented in [20]. The NP205 profile, which is part of the blade0 s structure, possesses a skeleton that is composed of quadratic polynomials and four-digit NACA profile thickness functions. The numerical results for different discharge values are used to determine the cavitation characteristics and sensitivity curves for the Francis pump-turbine.

19

! because v and φ are not dependent on θ and ! ! ! v ¼  ð i θ =RÞ  ∇Ψ , ∇  v ¼ 0 and ∂Ψ =∂θ ¼ 0. ! The components of velocity v are given by the following equations: vZ ¼

∂φ ; ∂Z

vZ ¼

1 ∂Ψ ; R ∂R

∂φ ∂R

vR ¼

vR ¼ 

ð3Þ 1 ∂Ψ R ∂Z

ð4Þ

Resolution of the problem with Dirichlet and Neumann boundary conditions for the stream function Ψ is achieved by integrating Eq. (2) using the FEM. To generalize the problem, it is useful to consider a dimensionless approach in the stream function, which uses the following variable interchange 1 Z n ¼ ZLax ;

1 Rn ¼ RLax

ð5Þ

and function

Ψ n ¼ 2π Q  1 Ψ

ð6Þ

where Lax is the axial extension of the analysis domain, Q is the discharge equal to 2π(Ψ  Ψ0) and the value of Ψ0 is zero. In this case, the Stokes equation for the stream function Ψn is expressed as ∂2 Ψ

n

∂Z n2

þ

∂2 Ψ

n

∂Rn2

1 ∂Ψ ¼0 Rn ∂Rn n



ð7Þ

and the dimensionless components of the velocity vn are given by vnZ n ¼

1 ∂Ψ ; Rn ∂Rn n

vnRn ¼ 

1 ∂Ψ Rn ∂Z n

n

ð8Þ

The equations that express the connection between the velocity dimensionless components and the dimensional components are expressed as follows:

2. Numerical model

vnZ n ¼ 2π L2ax Q  1 vZ ;

A two-dimensional FEM-DRM, which simulates flow velocity and pressure distributions on the blade of a Francis-type reversible radial-axial hydraulic machine runner in relative rotational motion, is presented in the next section.

Because the problem requires mixed limit conditions, Fig. 1 shows the geometric shape of the analysis domain and the limited conditions on its boundary for both the dimensional approach and the dimensionless approach. Function Ψn is globally approximated on Ωn as

2.1. Resolution of fluid axial-symmetric motion with the FEM

Ψ n ¼ anα Ψ nα ;

A method for solving the potential axial-symmetric motion of an inviscid incompressible fluid through the a runner of reversible hydraulic machinery without a blade, which operates as a pump, is developed. Because the fluid motion is axial-symmetric, the use of a cylindrical coordinate system, such as (R, θ, Z), is required. Previous assumptions indicate the potential of axial-symmetric ! ! motion, which suggests ∇  v ¼ 0, where v ¼ ∇φ. The potential φ of velocity ! v is stationary (∂φ/∂t¼0) because the absolute ! motion is considered stationary ð∂ v =∂t ¼ 0Þ. ! Considering that the fluid is incompressible ð∇ U v ¼ 0Þ, the following Stokes equations are obtained – with the velocity potential ∂2 φ

φ

∂2 φ

1 ∂φ ¼0 þ þ ∂Z 2 ∂R2 R ∂R

– with the stream function ∂2 Ψ ∂Z

2

þ

∂2 Ψ ∂R2



1 ∂Ψ ¼0 R ∂R

ð1Þ

Ψ ð2Þ

vnRn ¼ 2π L2ax Q  1 vR

α ¼ 1; :::; NG

ð9Þ

ð10Þ

where NG is the number of global nodes on domain Ωn, which resulted from its discretization of the finite elements, anα are the n n global interpolation functions on Ωn; and Ψ α is the value of Ψ in the global node α. The Galerkin method is applied as follows [15,21]: ! Z n n n ∂2 Ψ ∂2 Ψ 1 ∂Ψ n þ  ð11Þ anα dΩ ¼ 0 n n n2 ∂Rn2 R ∂R Ωn ∂Z The domain Ωn is axial-symmetric; thus, dΩn ¼RndθdZndRn. Using integration by parts in Eq. (11), a system of linear equations is obtained in its global form Dnαβ Ψ β ¼ F nα ; n

α; β ¼ 1; :::; NG

where coefficients Dnαβ are expressed as ! Z Z n n ∂anβ ∂anα ∂aβ ∂anα ∂aβ n n n1 n Dnαβ ¼ þ Ω þ 2 R a d α n dΩ n n ∂Rn ∂Rn ∂R Ωn ∂Z ∂Z Ωn and the free terms are expressed as ! Z n n ∂Ψ n ∂Ψ n n n n F nα ¼ ann n þ n α dΓ n Z ∂Rn R Γ n ∂Z

ð12Þ

ð13Þ

ð14Þ

20

A. Iosif, I. Sarbu / Engineering Analysis with Boundary Elements 41 (2014) 18–27 ne where the coefficients DNM and the free terms F nNe are computed with the following equations:  Z  ne ne Z ne ∂aN ∂aM ∂aN ∂ane ∂ane ne ne ne ¼ þ n Mn dΩ þ 2 Rn  1 anN Mn dΩ DNM n n ∂R ∂R ∂R Ωne ∂Z ∂Z Ωne ð17Þ

ne ¼ FN

Z Γ ne

! ne ne ∂Ψ ∂Ψ ne n n e n þ n ann n n N dΓ ∂Z n Z ∂Rn R

ð18Þ

e are numerically evaluated using a Gauss The coefficients DnNM quadrature formula [21] n

n

ne DNM ¼ 32  1 ∑ ∑ wi wj f NM ðζ i ; ηj Þ

ð19Þ

i¼1j¼1

where i is the number of Gauss integration points from the internal finite element, ζi and ηj are the natural coordinates of the integration point, and wi and wj are the weight factors. Considering the limit conditions for the stream function, the free terms F nNe are zero on Γne and Γn. Knowing the local values of ne ne Ψ Nne , DNM , and F N enables the computation of global values using Boolean matrices: E

Ψ nα ¼ ∑ ΔeNα ΨNne e¼1 E

ne Dnαβ ¼ ∑ DNM ΔN α ΔM β

e

e

e¼1 E

F nα ¼ ∑ F nNe ΔNα e

ð20Þ

e¼1

where ΔNα , ΔMβ are the elements of the Boolean matrix Δ , which has the dimensions NL  NG. The components of velocity vne on the finite element is computed in the gravity center of the finite element (for ζ ¼ 0, η ¼ 0) using Eqs. (15) and (21): e

e

e

vnZen ¼ 4α0 1 a2 1 AN2 Ψ N

ne

vRnen ¼ 4α0 1 a2 1 AN1 Ψ N

ne

ð21Þ

The velocity0 s value vn in the nodes is computed using the arithmetical average of the components that pertain to the neighboring finite elements: vn ¼ ðvnZ2n þ vRn2n Þ1=2

ð22Þ n

Fig. 1. Analysis domain and limit conditions for (a) the dimensional case and (b) the dimensionless case.

where anβ are the global interpolation functions on Ωn; Ψ β is the value of Ψn in the global node β; ann α are the global interpolation functions on Γn; and nnZ n , nnRn are the cosines of the angles between the normal nn to Γn boundary and the axis OZn, respectively the axis ORn. Establishing the global system of equations implies knowledge of the finite element equation. The analysis domain has been discretized into a number E of isoparametric linear finite elements Ωne with boundary Γne. The function Ψn is locally approximated on each finite element as n

Ψ ne ¼ anN Ψ nNe ;

N ¼ 1; :::; NL

ð15Þ

ne where aN are the local interpolation functions, Ψ N is the value of Ψne in the local node N, and NL¼4 is the number of local nodes on Ωne. If the Galerkin method is applied to a finite element, the form of the local finite element equation is ne

DNM Ψ M ¼ F N ; ne

ne

ne

N; M ¼ 1; :::; NL

ð16Þ

Knowing the components of the velocity v in the global nodes, is n possible to compute the values of ϕ in these nodes using the following equation Z φn ¼ φn0 þ vnZn dZ n þ vnRn dRn ð23Þ c

And imposing ϕ0 ¼ 0 on BC. The authors have developed a computer program in the FORTRAN programming language for PC compatible microsystems. Using this program, the hydrodynamic field (Fig. 2) defined by the streamlines Ψn ¼const and the equipotential lines φn ¼const are obtained. The velocity and pressure distributions (Figs. 3 and 4) along the streamlines are obtained considering the following data: Lax ¼2.027 m, an1 ¼ 0:666, vAB ¼ 12.61 m/s. The analysis domain was discretized for 492 finite elements and 546 global nodes. The velocity vn , which corresponds to a point on the streamline, represents the ratio between two dimensionless velocities n

vn ¼

vn vnAB

ð24Þ

where vnAB ¼ 2ðan1 Þ  2 ¼ 4:5

ð25Þ

A. Iosif, I. Sarbu / Engineering Analysis with Boundary Elements 41 (2014) 18–27

21

Fig. 4. Pressure distribution along the streamlines.

associated with the AB boundary and the second point is the current boundary. The point for which the maximum value for vn and the minimum value for p is obtained is identified in Figs. 3 and 4. The stream function Ψn ¼1 and the value sn ¼ 0:46 of the streamline curve arch are vnmax ¼ 1:364 and p ¼  0:86. These results indicate a higher probability of the cavitation phenomenon. For solving the rotational motion on the stream surface, only the streamlines from the hydrodynamic field in the meridian plane are retained (Fig. 5). The entry and exit edges of the blade are plotted in this plane; the latter is parallel to the axis OZn and is located at radius RnP ¼ 1:233.

Fig. 2. Hydrodynamic field.

2.2. Flow analysis around radial–axial profile cascades with the DRM 2.2.1. Basis of the DRM The objective of the DRM proposed by Brebbia and Nardini [22,23] and developed by Novak, Partridge and Wrobel [18,24,25] is based on the BEM: the transformation of domain integrals into boundary integrals. The elements that constitute the basis of the method are expressed as ∇2 u ¼  b

ð27Þ

where b is a dependent function of (x, y, u, t) or of (x, y, u) when time t is missing. Solution u of Eq. (27) can be obtained from solution u of the Laplace equation and a particular solution u^ as follows: u ¼ u þ u^

ð28Þ

In Refs [18,26], the use of a limited series particular solutions and Eq. (29) for function b is proposed NþL

b ¼ ∑ F j αj

Fig. 3. Velocity distribution along the streamlines.

j¼1

and an1 is the dimensionless radius that corresponds to the entry in the considered runner (Fig. 1b). The dimensionless velocity v for every streamline point is obtained by multiplying vn by vAB. To obtain the pressure distribution along the streamlines, (Fig. 4) the following equation is employed: p ¼ 1  vn2

ð26Þ

Eq. (26) is obtained using Bernoulli0 s theorem for the potential stationary motion for two streamline points: the first point is

ð29Þ

where Fj is the function associated with point j and αj are the α vector elements. Because the boundary contains N nodes and the Ω domain contains L internal points (Fig. 6), N þL values are obtained for Fj and αj. Beginning with the basic equation NþL

∇2 u ¼  ∑ ð∇2 u^ j Þαj j¼1

ð30Þ

to which the standard BEM technique is applied to Γ boundary0 s discretization of the linear elements, the following matrix equation is

22

A. Iosif, I. Sarbu / Engineering Analysis with Boundary Elements 41 (2014) 18–27

Fig. 7. Radial-axial cascade on the stream surface.

The computation of elements Hij and Gij, which correspond to matrices H and G, is equivalent to the computation of linear boundary elements. After implementation of the limited conditions, the matrix Eq. (31) produces a linear system of N þL equations with N þL unknowns, which consist of the values of function u and its normal derivative q in the nodes in which they are unknown. 2.2.2. Equation of relative rotational motion on a stream surface The fluid flow in a Francis-type runner is considered fluid and incompressible. The crossing of the stream surface with the runner blades, when operating as a pump, will determine a radial-axial cascade with profiles, as illustrated in Fig. 7. If the potential motion of ideal fluid, the relation between absolute velocity v, relative velocity w and peripheral velocity u, and the Oq1q2q3 orthogonal curve-linear coordinate system are considered, the stream surface motion equation [20] is obtained: ∂wð1Þ ∂ ½RR0 1 wð2Þ   þR0 1 Rf P ¼ 0 ∂q1 ∂q2 Fig. 5. Streamlines and entry and exit edges for the runner blades.

ð33Þ

where f P ¼ 2 ωP

"  2 #  ð1=2Þ dR dR 1þ dz dz

ð34Þ

where w(1) and w(2) are the relative flow velocity components of motion on the stream surface, ωP is the angular velocity; R0 is the radius corresponding to the origin of curvilinear coordinates system; and z, R is the axial coordinate and radius variables in the case of axisymmetric motion. Fig. 6. Boundary and internal points from the Ω domain.

2.2.3. Conformal mapping of the radial-axial cascade Because Eq. (33) is not required in numerical calculations, the following change of variables (x, y) is necessary: Z q1 x¼ R0 R  1 dq1 ; y ¼ q2 ð35Þ

obtained ^ α Hu  Gq ¼ ðGQ^ HUÞ

ð31Þ

where u, q are the solution of its normal derivative vectors with size ^ have the size (Nþ L)  (NþL). To NþL, and matrices H, G, Q^ , H ^ a limited series that consists of compute the matrix elements Q^ , H, the particular solution u^ and its normal derivative q^ is proposed, in which the distance function r intervenes: u^ ¼

Eq. (35) realizes the geometric transformation of a stream surface in the associated plane and the transformation of the radial-axial cascade into a linear cascade. In numerical computations, these transformations are accomplished in the associated plane in which the linear cascade is located, as shown in Fig. 8; the following dimensionless variables are employed xn ¼

r2 r3 rm þ 2 þ þ :::: þ 4 9 ðm þ 2Þ2

   ∂x ∂y 1 r rm þ þ :::þ q^ ¼  x þ y ∂n ∂n 2 3 mþ2

0

ð32Þ

xn ; Ln

where Z xn ¼

yn ¼

sn

sn0

yn Ln

r n0 r n  1 dsn ;

ð36Þ

yn ¼ r n0 θ

ð37Þ

A. Iosif, I. Sarbu / Engineering Analysis with Boundary Elements 41 (2014) 18–27

in which

α¼F

1

23

α, λ are vectors, where

a;

λ ¼ F  1c

ð45Þ

In this equation a, c are vectors, and F–1 is the reciprocal of matrix F. By solving matrix Eq. (44) with the boundary conditions, Eq. (41) is performed through an iterative process, beginning with the Laplace equation solution. The results from the plane are transposed on the stream surface using the following relations wð pð

R0 ; Lax

rn ¼

R Lax

ð38Þ

in which sn0 is the dimensionless arch that corresponds to the origin O (Fig. 7), Lax is the dimensional axial extension of the analysis domain, and Ln is the cascade chord length in the associated plane for the dimensionless problem in the stream function Ψn solved in axisymmetric motion. 2.2.4. Solving the dimensionless differential equation of the stream function in the associated plane Considering Eq. (33) in the associated plane, the following differential dimensionless equation for the ψn stream function have been deduced ∂2 ψ n ∂xn2

þ

∂2 ψ n ∂yn2

1 dhðxn Þ ∂ψ n n  Þ  hðxn Þðr 0n  1 r n Þ2 f P ¼ 0 hðxn dxn ∂xn

ð39Þ

where f P ¼ 2ωnP n

"  n 2 #  1=2 dr n dr 1 þ dzn dzn

ð40Þ

and hðxn Þ is the thickness function of the variable fluid layer. Differential Eq. (39) is solved with the DRM considering the analysis domain shown in Fig. 8 and the following boundary conditions

ψ n ¼ 0 on JI

ð41Þ

n

where t 0 is the cascade step from the associated plane, nn is the external normal to the boundary, and β is the angle between the relative and peripheral velocities. Eq. (39) can be expressed as ∇2 ψ n ¼ b

ð42Þ

where b is the sum of the following two terms: n

a ¼ hðxn Þðr 0n  1 r n Þ2 f P ;

c¼

0

h ðxn Þ ∂ψ n Þ hðxn ∂xn

ð43Þ

If the DRM obtained for Eq. (39) is considered, the following matrix equation is obtained ^ GQ^ Þðα þ λÞ Hψn  Gq ¼ ðHU

Þ

¼ 1 ðR0 R  1 Þ2 wn2 þ ωP 2 R20  ðvm0 wnAM Þ  2 ½ðRR0 1 Þ2  1

ð47Þ

wn ¼

wn wnAM

ð48Þ

where wn is the dimensionless velocity computed in the associated plane, wnAM is the dimensionless upstream velocity, and vm0 is the dimensional meridian velocity that corresponds to the “O” origin of the curvilinear coordinate system. A computer program has been elaborated based on the developed numerical model. It was realized in the FORTRAN programming language for PC compatible microsystems.

3. Application of the numerical simulation model to a pump-turbine The proposed computational model is applied to reversible radial-axial hydraulic machinery, which operates as a pump, with the following fundamental characteristics: Hc ¼317 m, Qc ¼72.35 m3/s, Pc ¼250 MW, ηc ¼0.9, nsp ¼130, nc ¼ 300 rpm, AM AV Ln ¼ 1.476, βAM ¼1591, βc ¼211, βAV ¼ 1551, β c ¼251, βsup ¼241, and βc ¼ 661. The number of runner blades is ZP ¼7 and the maximum relative thickness of the profile is d/l ¼0.5. The thickness function of the variable fluid layer is expressed as follows: hðxn Þ ¼ 57:74xn5  81:38xn4 þ 38:67xn3  7:27xn2  0:17xn þ 0:998 ð49Þ

ψ n ¼ t n0 on O0 K ∂ψ n AM AM ¼  ctg β c on AB n ¼ ctg β

∂n ∂ψ n AV AV ¼  ctg β ¼ ctg βc on HG ∂nn ψ n ¼ ψ n AJ; IH þ t n0 on BO0 ; KG  n ∂ψ n ∂ψ ¼  on BO0 ; KG ∂nn ∂nn AJ; IH

ð46Þ

and by obtaining the velocity and pressure distributions (wðÞ , pðÞ ) on the profile from the radial-axial cascade of the reversible runner blades. The dimensionless velocity wn is expressed as

Fig. 8. Analysis domain from the associated plane.

r n0 ¼

¼ R0 R  1 w n

Þ

ð44Þ

To solve the matrix in Eq. (44), the analysis domain is discretized into linear elements on the boundary and the internal points are established. The numerical computation is performed for N ¼ 43 nodes of the boundary and L ¼20 internal points (Fig. 9). To iteratively solve Eq. (44) using the DRM with a precision of 2  10  4, 18 iterations are required. The flow velocity and pressure distributions on a reversible runner blade are simulated with a computer program for the discharge Qc ¼ 72.35 m3/s and the medium stream surface generated by the streamline Ψn ¼0.6, as illustrated in Figs. 10 and 11. This representation was performed along the LnLOX dimensionless loxodrome. For a comparison of the results, Figs. 10 and 11 also present the velocity and pressure distributions computed with the FEM. Fig. 12 shows computational grids for the linear cascade of xn o yn in the associated plane. The analysis domain is discretized for 209 linear isoparametrical finite elements and 240 global nodes. From the representation of velocity and pressure distributions (Figs. 10 and 11), the following values on profile intrados (soffit) were obtained using the DRM with wð Þ ¼1.56, pð Þ ¼ –1.17 and the FEM with wð Þ ¼1.46 pð Þ ¼  1. Satisfactory values were obtained with the two numerical methods.

24

A. Iosif, I. Sarbu / Engineering Analysis with Boundary Elements 41 (2014) 18–27

Fig. 11. Pressure distribution on the runner blade.

Fig. 9. Discretization of the analysis domain for the DRM.

Fig. 12. Discretization of the analysis domain for the FEM.

Fig. 10. Velocity distribution on the runner blade.

4. Cavitation characteristics and sensitivity curves The reversible hydraulic machine operates at random discharge values of Qx, which are either different than or equivalent to the

computed discharge. Considering the unfavorable behavior of the machine during the pump regime compared with the turbine regime [1,27,28], it is necessary to determine the cavitation characteristics sPx ¼ f(Qx/Qc) and the cavitation sensitivity curves kp max, x ¼f(Qx/Qc) during the runner design stage. These theoretical characteristics are determined based on the coefficient of minimal pressure distributions obtained for Qx ¼ [0.8, 1.0, 1.2, 1.4]  Qc and Ψn ¼0.6. Fig. 13 presents the Francis

A. Iosif, I. Sarbu / Engineering Analysis with Boundary Elements 41 (2014) 18–27

25

Fig. 13. Hydraulic circulation of the reversible machine, which operates as a pump. Fig. 16. Pressure distribution on the hydraulic machine0 s blade.

Fig. 14. Velocity triangles that correspond to point 0.

Fig. 17. Cavitation characteristics sPx ¼ f(Qx/Qc).

The pumping head HPx is computed as H Px ¼

Fig. 15. Velocity distribution on the hydraulic machine0 s blade.

pump-turbine hydraulic circulation while operating as a pump, and Fig. 14 shows the velocity triangles that correspond to points 0 on the streamlines at the discharge Qx. In the case of Qx random discharge regimes, the velocity (Fig. 15) and pressure (Fig. 16) distributions are obtained using the DRM for one of the blades of the reversible hydraulic machine. To compute the pump cavitation coefficient sPx for the discharge Qx, the following equation is employed

sPx ¼ kp max ;x

w20x u2 v2 ∑hp0  Mx aMx DS  kux 0x þ 0x þ þ 2gH Px 2gH Px 2gH Px H PX H Px

where hydraulic losses ∑hp 0 incompressible fluid.

 Mx

ð50Þ

are zero for an inviscid

ηhP

gð1 þ pÞ

 u2P u2P 

Qx ctg β2P π D2P b2P

 ð51Þ

where p is the influence coefficient for a finite number of blades ZP, (ZP ¼ 7, p ¼0.3836) and ηhP ¼ 0.96 is the hydraulic efficiency of the runner. In Eq. (50), the factor kp max, x is equivalent to the value  pð Þmin ;x that was obtained from the pressure distributions for discharge Qx, and the absolute velocity (v0x), relative velocity (w0x) and peripheral velocity (u0x) that correspond to point 0 in Fig. 13, with respect to origin O of the curvilinear coordinate system Oq1q2q3 in Fig. 7. The numerical results obtained with the DRM illustrate the cavitation characteristics (Fig. 17) and of the cavitation sensitivity curves (Fig. 18) with the respective characteristics determined from the FEM. In Fig. 17, the value of sP of 0.145, which was obtained with the statistical equation Qx ¼Qc recommended by Siervo [29], was also established. This value is closer to the sP values computed with the DRM and the FEM for the intrados (soffit). Figs. 17 and 18 show that an increase in discharge Qx causes an increase in the values of sPx and kp max, x, which produces an unfavorable runner from a cavitational point of view. These operating conditions enable the runner to be developed inside the cavitation and reduce the energy efficiency.

26

A. Iosif, I. Sarbu / Engineering Analysis with Boundary Elements 41 (2014) 18–27

Table 1 The coordinates of point M0 s location.

Fig. 18. Cavitation sensitivity curves kp

Method

xnM

xnM

snM

Z nM

RnM

ku

sP

DRM FEM

0.220 0.238

0.325 0.351

0.498 0.526

0.494 0.520

0.577 0.587

0.23 0.27

0.216 0.175

max, x ¼ f(Qx/Qc).

A computational method for the cavitation coefficient sPx and the M point location, in which the pressure is expressed by pð Þmin ;x for Ψn ¼0.6 and Qx ¼Qc ¼ 72.35 m3/s, is presented in the next section. For the case in which Qx/Qc ¼1, Eq. (50) can be formulated as

sP ¼ kp max

w20 u2 v2 ∑hp0  M aM DS  ku 0 þ 0 þ þ 2gH P 2gH P 2gH P HP HP

ð52Þ

where the peripheral velocity coefficient is given by  2  n 2 RM RM ku ¼ 1 ¼ 1 R0 Rn0

ð53Þ

and aM DS ¼ ð0:8625  Z nM ÞLax n

ð54Þ n

where radius RM of point M and radius R0 of point 0 are dimensionless and RM ¼ RnM Lax and R0 ¼ Rn0 Lax are dimensional. In Eq. (52), the hydraulic losses ∑hp 0  M are zero and the values for the velocities are v0 ¼12.9 m/s, u0 ¼33.45 m/s, and w0 ¼35.3 m/s. To establish the location of point M, the following approximation polynomials have been determined: 1 ¼ 0:10827sn6 1:5184sn5 þ 4:5641sn4  Rn  3:809sn3 þ 0:27393sn2 þ 0:10227sn þ 1:9152

ð55Þ

Z n ¼ 0:023713sn5 þ 0:27655sn4  1:1190sn3 þ þ 0:80934sn2 þ 0:80887sn þ 0:011191

ð56Þ

sn ¼ 0:53727xn4 þ 0:029893xn3 þ 0:02595xn2 þ þ 1:0012xn þ 0:1624

ð57Þ

For sn ¼ 0, the result is Rn0 ¼ 0:522 and the dimensionless abscissa of point M is given by the relation xnM ¼ Ln xnM n

ð58Þ n

n

where xM is the abscissa of point M in plane y ox . The computation results are summarized in Table 1. Fig. 19 illustrates the location of point M that is obtained using the numerical DRM and the FEM.

Fig. 19. Location of point M.

Minor differences between the results obtained with the DRM and the results obtained with the FEM are observed, which can be attributed to the computational errors generated by these methods. The value of the pump cavitation coefficient sP obtained using the proposed model is comparable with the 3D simulation result performed by other researchers for similar input data using FLUENT commercial software. Stuparu et al. [5] analyzed the cavitational behavior of a large storage pump with two identical stages and considered the following fundamental characteristics: Qc ¼1 m3/s, Hc ¼159, 5 m, and ZP ¼ 5. Using the FLUENT program, the values for the incipience cavitation coefficient si, which are presented in Table 2, were determined. Because bubbles are initially formed at incipience cavitation, only si ¼0.44 for Qc ¼1.02 m3/sE1 m3/s is of interest. The result sinst ¼0.23 is obtained from the graphical illustration (Fig. 7 [5]) of the variation in the pumping head (H), which is dependent on the installation cavitation coefficient (sinst) for Qc ¼QOP3 ¼1.02 m3/s and Hc ¼159.5 m.

A. Iosif, I. Sarbu / Engineering Analysis with Boundary Elements 41 (2014) 18–27

Table 2 Values of the incipience cavitation coefficient for the operating points. Operating points

Qx [m3/s]

Qx/Qc [dimensionless]

si [dimensionless]

OP1 OP2 OP3 OP4 OP5

0.8 0.9 1.02 1.1 1.2

0.8 0.9 1.0 1.1 1.2

0.32 0.38 0.44 0.49 0.53

The pump operates in acoustic detectable cavitation when

sP ¼ sinst ¼ 0.23. Therefore, this value corresponds with the values

obtained using the proposed model, which is based on the DRM (sP ¼ 0.216) and the FEM (sP ¼0.175). Point M of minimal pressure is located in the vicinity of the area where the liquid flow on the hydraulic machine blade exists in a pumping state, which has been experimentally identified by Anton [27]. 5. Conclusions The authors developed a methodology for solving the 3D problem of flow in reversible hydraulic machinery by transforming it into two 2D problems to determine velocity and pressure distributions on a pump-turbine runner blade using FEM-DRM coupling, in which an inviscid fluid is assumed. For the numerical simulation of velocity and pressure distributions, it is necessary to solve a mixed-boundary conditions problem in the associated plane. A numerical computational model, which is based on the DRM, has been developed to solve the partial derivative in Eq. (39) by imposing Dirichlet and Neumann boundary conditions on the analysis field boundary in the associated plane. This model comprises an innovative approach to solving the problem. The modality of solving the problem with the DRM can serve as a starting point in the case of a viscous fluid. A theoretical study of the cavitation characteristics and sensitivity curves and the establishment of the minimal pressure point location yields important results for design engineers in this field. References [1] Coutier-Delgosha O, Fortes-Patella R, Rebound JL, Hofmann M, Stoffel B. Experimental and numerical studies in a centrifugal pump wit two-dimensional curved blades in cavitating conditions. J Fluids Eng 2003;125:970–8. [2] Cudina M. Detection of cavitation phenomenon in a centrifugal pump using audible sound. Mech Syst Signal Process 2003;17(6):1335–47. [3] Esxaler X, Egusquiza E, Farhat M, Avellan F, Coussirat M. Detection of cavitation in hydraulic turbines. Mech Syst Signal Process 2006;20:983–1007. [4] Pouffary B, Patella RF, Rebound JL, Lambert PA. Numerical simulation of 3D cavitating flows: analysis of cavitation head drop in turbomachinery. J Fluids Eng 2008;130:061301-1–10.

27

[5] Stuparu A, Susan-Resiga R, Anton LE, Muntean S. A new approach in numerical assessment of the cavitation behavior of centrifugal pumps. Int J Fluid Mach Syst 2011;4(1):104–13. [6] Anton I. Hydraulic turbines. Publishing House Facla. 1979 (in Romanian). [7] Dring RP, Joslyn HD, Hardin LW, Wagner JH. Turbine rotor–stator interaction. ASME J Eng Power 1982;104:729–42. [8] Kolšek T, Duhovnik J, Bergant A. Simulation of unsteady flow and runner rotation during shutdown of an axial water turbine. J Hydraul Res 2006;44(1):129–37. [9] B.S.K. Naidu, Computerized design development and model testing of a high head Francis turbine runner. In: Proceedings of the 15th IAHR symposium. Belgrade, Yugoslavia, Q4; (1990) p. 1–11. [10] Nova RA, Hearsey RM. A nearly three-dimensional interblade computing system for turbo-machinery. ASME J Fluids Eng 1994;99(3):154–66. [11] Qian Y, Suzuki R, Arakawa C. Analysis of the inlet reverse flows in a pump turbine using 3-D viscous numerical techniques. Hydraulic machinery and cavitation, 1. Netherlands: Kluwer Academic Publishers; 1996; 230–8. [12] Song CS, Chaugsi C, Ikohagi T, Sato J, Shinmei K, Tani K. Simulation of flow through pump–turbine. Hydraulic machinery and cavitation, 1. Netherlands: Kluwer Academic Publishers; 1996; 277–84. [13] Abdallah S, Henderson RE. Improved approach to the streamline curvature method in turbo-machinery. ASME J Fluids Eng 1987;109(3):213–7. [14] S.F. Timouchev, K.P. Illichov, 2D numerical simulation of unsteady pressure field in centrifugal pumps. In: Proceedings of International INSE symposium. Lyon, France, 3; (1995). [15] Reddy JN. An Introduction to the finite element method. New York: McGraw-Hill; 1993. [16] Katsikadelis J. Boundary element method, theory and application. Oxford, UK: Elsevier; 2002. [17] Sarbu I. Numerical modeling and optimizations in building services. Timisoara: Politehnica Publishing House; 2010 (in Romanian). [18] Partridge WP, Brebbia AC, Wrobel CL. The dual reciprocity boundary element method. Southampton: Computational Mechanics Publications; 1992. [19] Carte IN. Determination of the incompressible fluid motion round radial–axial profile cascades through the finite element method. Appl Mech 1987;32(6): 639–49. [20] A. Iosif, Radial–axial reversible hydraulic machinery, Ph.D. thesis, Politehnica University of Timisoara, Romania, 1998. [21] I.N. Carte, Contributions to the study of radial–axial cascades and applications for Francis turbines design, Ph.D. thesis, Polytechnic Institute of Timisoara, Romania, 1986. [22] Brebbia AC, Nardini D. Dynamic analysis in solid mechanics by an alternative boundary element procedure. Int J Soil Dyn Earthq 1983;2(4):228–33. [23] D. Nardini, A.C. Brebbia, Boundary integral formulation of mass matrices for dynamic analysis, In: Boundary element research, Springer-Verlag, Berlin, 2 (1985) 191-208. [24] Novak AJ, Brebbia AC. The multiple reciprocity method: a new approach for transforming BEM domain integrals to the boundary. Eng Anal 1989;6(3): 164–7. [25] Partridge WP, Brebbia AC. Computer implementation of the BEM dual reciprocity method for the solution of general field equations. Commun Appl Numer Methods 1990;6(2):83–92. [26] Wrobel CL, Brebbia AC. The dual reciprocity boundary element formulation for nonlinear diffusion problems. Comput Methods Appl Mech Eng 1987;65(2): 147–61. [27] I. Anton, Cavitation, Romanian Academy Publishing House, Bucharest, 1985 (in Romanian). [28] Nouri NM, Eslamdoost A, Shienejad A, Moghimi M. Calculation of hydrodynamic drag force exerted on a supercavitating slender body using boundary element method. WSEAS Trans Fluid Mech 2006;1(6):566–71. [29] Siervo F, Lugaresi A. Modern trends in selecting and designing reversible Francis pump–turbines. Water Power Dam Constr 1980;5(1980):234–45.