Mentzos M Etal_2005_cfd Predictions Of Flow Through A Centrifugal Pump Impeller

1st International Conference on Experiments/Process/System Modelling/Simulation/Optimization 1st IC-EpsMsO Athens, 6-9 July, 2005 © IC-EpsMsO

CFD PREDICTIONS OF FLOW THROUGH A CENTRIFUGAL PUMP IMPELLER 1

1

1

Michalis D. Mentzos , Andronicos E. Filios A.E.*, Dionisios P. Margaris and Dimitrios G. Papanikas 1 Fluid Mechanics Laboratory, Mechanical Engineering and Aeronautics Department University of Patras GR-26500 Patras, Greece e-mail: [email protected] *

Fluid Mechanics & Turbomachines Laboratory, Mechanical Engineering Department, School of Pedagogical and Technological Education, GR-141 21 Athens, Greece e-mail: [email protected]

Keywords: Centrifugal pump, Impeller, Computational Fluid Dynamics. Abstract. With the aid of computational fluid dynamics, the complex internal flows in water pump impellers can be well predicted, thus facilitating the design of pumps. This paper describes the simulation of the flow into the impeller of an experimental centrifugal pump. The commercial three-dimensional Navier-Stokes code Fluent, with a standard k-ε two-equation turbulence model is used to simulate the problem under examination. In the calculation, the finite-volume method along with a structured grid system is used for the solution procedure of the discretized governing equations for this problem. The calculation predicts the flow pattern and the pressure distribution in the untwisted blade passages and finally the overall performances and the head-capacity curve are discussed. 1 INTRODUCTION Pump designers are continually challenged to provide machines that operate more efficiently, quietly, and reliably at lower cost. Key to design a hydraulic turbomachine is a detailed understanding of the internal flow within its stationary and rotating passages and therefore the calculability of its performance during design and off-design conditions. With the aid of the Computational Fluid Dynamics (CFD) approach, the complex internal flows in water pump centrifugal impellers, which are not fully understood yet, can be well predicted and therefore establishing the CFD as a key tool for pump designers. The use of CFD tools in turbomachinery industry is quite common today since many tasks can numerically be solved much faster and cheaper than by means of experiments. Commercial software Fluent®, CFX-Tascflow®, StarCD® and Fine/Turbo®, is increasingly using to study pump design and off-design performance. Various researchers have considerably contributed in revealing the flow mechanisms inside centrifugal impellers aiming to the design of high performance centrifugal turbomachines. The reported works by Eckard[1], Johnson and Moore[2], Kjork and Lofdahl[3], Denton[4], Dawes[5], Casey et al[6], Bansod and Rhie[7], Krain and Hoffman[8], Farge and Johnson[9] and Zhang et al[10] are an indicative collection of research efforts on the computation and the experimental verification of the flowfields within centrifugal impellers. On the numerical simulation of the flowfield of a centrifugal impeller, several algorithms have been proposed and developed, but there is still no approach that is fully robust regarding the numerical and modeling accuracy as well as efficiency. Pressure-based methods initially developed for the incompressible flow regime, obtain the pressure field via a pressure or a pressure correction equation which is formulated by manipulating the continuity and momentum equations. The solution procedure is conventionally sequential in nature, and hence can more easily accommodate a varying number of equations depending on the physics of the problem involved, without the necessity of reformulating the entire algorithm. Lakshminarayana[11], Rodi et al[12], Thakur et al[13], provides a review of the techniques that are useful as an assessment of the state of the art. The present paper summarizes the under progress research work on the prediction of the flowfield within an experimental centrifugal pump as well as its performance curves. Towards to study the design and off-design performance of the pump, flow and pressure fields are analysed numerically and experimentally for reasons of direct comparison. The experimental pump has a shrouded impeller of constant width rotating inside a cylindrical casing. The impeller have nine untwisted blades backward facing with a circular arc camber line and its cad model along with the main geometry characteristics are shown in figure 1. In the lack of completing the experimental study, only the computational study of the flow through the impeller with the use of commercial

Michalis D. Mentzos, Andronicos E. Filios A.E., Dionisios P. Margaris and Dimitrios G. Papanikas ®

software Fluent is discussed. The software package Fluent® can predict laminar flow, turbulent flow, and heat transfer. It has been widely used in the field of turbomachinery, and the simulation results have been proven by Sun and Tsukamoto[14] and Gonzalez et al.[15] to be reliable. Fluent® overcomes the meshing difficulties that arise in complex geometry by using a powerful CAD-based preprocessor, Gambit®. The authors [16] have already used Fluent® to study three-dimensional turbulent flow through a commercial water-pump during design and off-design conditions. The flow and pressure field through the impeller results from the solution of the fully three-dimensional incompressible Navier-Stokes equations, including the centrifugal force source. Turbulence is simulated with the standard RNG k-ε model. Although grid size is not adequate to investigate local boundary layer variables, global ones are well captured. For such calculations, wall functions, based on the logarithmic law, are used. The pressure-velocity coupling is calculated through the SIMPLE algorithm. Second order, upwind discretizations is used for convection terms and central difference schemes for diffusion term.

Denomination Impeller diameters Blade angles Impeller widths Number of blades Blade thickness Blade curvature radius

Value D1= 62 mm, D2 = 190 mm β1 = 26o, β2 = 49o b1 = 9 mm, b2 = 9 mm z=9 t = 5 mm Rb= 110 mm

Figure 1: CAD-model of the centrifugal impeller under study.

2 MATHEMATICAL MODELING 2.1 Basic Equations For three-dimensional incompressible, unsteady flow, the continuity and momentum equations can be written in the rotating coordinate system as follows:

∂ρ + ∇ ⋅ (ρU ) = 0 ∂t

(1)

(

(

))

∂ρU + ∇ ⋅ (ρU ⊗ U ) = ∇ ⋅ − pδ + µeff ∇U + (∇U )T + S M ∂t

(2)

where vector notation has been used, ⊗ is a vector cross product; U is the velocity; p is the pressure; ρ is the density; δ is the identity matrix; SM is the source term and µeff is the effective viscosity coefficient. For flows in a rotating frame of reference at a constant rotational speed Ω, the Coriolis effect is taking into account, too. In this case, with r the location vector, S M = − ρ [2 Ω ⊗ U + Ω ⊗ (Ω ⊗ r )]

(3)

2.2 RNG k-ε turbulence model The RNG-based k-ε turbulence model is derived from the instantaneous Navier-Stokes equations, using the

Michalis D. Mentzos, Andronicos E. Filios A.E., Dionisios P. Margaris and Dimitrios G. Papanikas

“renormalization group” (RNG) technique. The analytical derivation results in a model with constants different from those in the standard k-ε model, and additional terms and functions in the transport equations for k and ε.

.

The RNG k-ε model has a similar form to the standard k-ε model The values of k and ε come directly from the differential transport equations for the turbulence kinetic energy and turbulence dissipation rate

∂ ∂ ∂ ∂k (ρ k ) + ( ρ kui ) = (a µ ) + G k − ρε ∂t ∂x i ∂ x j k eff ∂ x j

(4)

∂ ∂ ∂ ∂ε ε ε2 ( ρε ) + ( ρε ui ) = ( aε µ eff ) + C 1ε G k − C 2 ε ρ − Rε k k ∂t ∂x i ∂x j ∂x j

(5)

and

where Gk is the production of turbulence kinetic energy due to the mean velocity gradients; αk, and αε are the inverse effective Prandtl numbers for k and ε, respectively; C1ε and C2ε are constants with C1ε=1,42 and C2ε=1,68, respectively. The term Gk, representing the production of turbulence kinetic energy, is modeled identically for the standard, RNG, and realizable k-ε models. From the exact equation for the transport of k, this term may be defined as

Gk = − ρ ui′u′j

∂u j

(6)

∂x i

The scale elimination procedure in RNG technique results in a differential equation for turbulent viscosity

⎛ ρ 2k d⎜ ⎜ εµ ⎝

⎞ vˆ dvˆ ⎟⎟ = 1, 72 3 (vˆ − 1 + C v ) ⎠

(7)

with Cv ≈100 and

νˆ =

µ eff

(7a)

µ

In the high - Reynolds numbers limit, equation (7) is integrated to obtain an accurate description of how the effective turbulent transport varies with the effective Reynolds number (or eddy scale), allowing the model to better handle low-Reynolds-number and near-wall flows

µt = C µ ρ

k2

ε

with Cµ =0,0845, derived using RNG theory. It is interesting to note that this value of Cµ is very close to the empirically-determined value of 0,09 used in the standard k-ε model. The inverse effective Prandtl numbers for k and ε, αk and αε respectively are computed by

α − 1,3929 α 0 − 1,3929

0 ,6321

α + 2,3929 ⋅ α 0 + 2,3929

0 ,3679

=

µ mol µ eff

(8)

where α0 =1,0. In the high Reynolds number limit (µmol/µeff<<1), αk=αε≈1,393 The main difference between the RNG and standard k-ε models lies in the additional term in the ε equation given by Rε =

C µ ρη 3 [1 − (η / η 0 )] ε 2 1 + βη 3

k

(9)

where η≡Sk/ε, η0=4,38, β=0,012 The effects of Rε-term in the RNG ε-equation can be seen more clearly by rearranging equation (5).

Michalis D. Mentzos, Andronicos E. Filios A.E., Dionisios P. Margaris and Dimitrios G. Papanikas

Using equation (9), the third and fourth terms on the right-hand side of equation (5) can be merged, and the resulting ε-equation can be rewritten as

∂ ∂ ∂ ∂ε ε ε2 ( ρε ) + ( ρε ui ) = ( aε µ eff ) + C 1ε G k − C 2*ε ρ ∂t ∂xi ∂x j ∂x j k k

(10)

where C2ε* is given by

C 2*ε

C µ ρη 3 ⎛⎜ 1 − η ⎞⎟ η0 ⎠ ⎝ = C 2ε + 3 1 + βη

(10a)

In regions where η<η0, the Rε-term makes a positive contribution, and C2ε* becomes larger than C2ε. In the logarithmic layer, for instance, it can be shown that η≈3, giving C2ε*≈2, which is close in magnitude to the value of C2ε=1,92 in the standard k-ε model. As a result, for weakly to moderately strained flows, the RNG model tends to give results largely comparable to the standard k-ε model. In regions of large strain rate (η>η0), however, the Rε-term makes a negative contribution, making the value of less C2ε* than C2ε. In comparison with the standard k-ε model, the smaller destruction of ε augments ε, reducing k and, eventually, the effective viscosity. As a result, in rapidly strained flows, the RNG model yields a lower turbulent viscosity than the standard k-ε model. Thus, the RNG model is more responsive to the effects of rapid strain and streamline curvature than the standard k-ε model, which explains the superior performance of the RNG model for certain classes of flows. Equations (1), (2), (4) and (5) form a closet set of nonlinear partial differential equations governing the fluid motion. 2.3 Log-law wall functions There are large gradients in the dependent variables near the wall. It is costly to fully resolve the solution in this near-wall region as the required number of nodes would be quite large. Thus a common approach known as “wall functions” is applied to model this region. In the wall-function approach (Launder and Spalding[17]), the near-wall tangential velocity is related to the wall shear stress by means of a logarithmic relation, which can be written as follows:

u+ =

where

u+ =

1

κ

( )

ln y + + C

ut , + ρ∆yut and τw y = uτ = µ uτ ρ

(13)

(13a)

τw is the wall shear stress, ut is the known velocity tangent to the wall at a distance of ∆y from the wall, κ is the von Karman constant for smooth walls, and κ and C are constants, depending on wall roughness. However, this form of the wall-function equations has the problem that it becomes singular at separation points where the near-wall velocity, ut, approaches zero. In the logarithmic region, the alternative velocity scale, u*, can be used instead of u+:

u* = c µ1 / 4 k

(14)

This scale has the useful property of not going to zero if ut goes to zero (and in turbulent flow, k is never completely zero). Based on this definition, the following explicit equation for the wall shear stress is obtained:

τ w = τ visc

y* u+

(15)

where τvisc= µut/∆y; y*= ρu*∆y/µ; and u+ as given from equation (13). The recommended practice is to locate near-wall nodes such that y* is in the range of 20 to 50 for smooth walls. In the near-wall region, an estimate of the dissipation consistent with the log-law can be presented as

ε=

c µ3 / 4 k 3 / 2 u∆ y

(16)

Michalis D. Mentzos, Andronicos E. Filios A.E., Dionisios P. Margaris and Dimitrios G. Papanikas

The dissipation at the first interior node is set equal to this value. The boundary nodal value for k is estimated via an extrapolation boundary condition. The near-wall production of turbulent kinetic energy is derived to be

pk =

2 τ visc p k∗ µ

(17)

where 2

pk∗

⎛ y ∗ ⎞ du + =⎜ +⎟ ⎜ u ⎟ dy ∗ ⎝ ⎠

(18)

3 COMPUTATIONAL MODELING 3.1 Meshing Modeling the impeller numerically, ideally all nine impeller passages should be included in order to simulate the true flow field and detect any asymmetry of the flow. However, most numerical simulations of the flow field in turbomachines take advantage of the geometrical symmetry of the impeller and simulate one impeller passage only, as it is shown in figure 2. The computational grid is of structured type and it is generated with the Fluent preprocessor Gambit. The whole domain consists of three subdomains or zones, in a way that the density and the quality of the cells in local flow field regions can be suitably controlled and handled depending on pressure gradients and velocities. The fist and third zones are stationary while the second zone that incorporates the blade is moving with the applied rotational speed, i.e. 1450 rpm or 1950rpm. The first zone represents the suction or inlet pipe and the third zone is the discharge or outlet portion where the flow is fully developed with a less possible reacting outlet boundary condition. Structured hexahedral cells are used in the whole domain with 41.211 cells for the inlet zone, 48.000 cells for the impeller 14.040 cells for the outlet zone. The size of the resulting cells is not adequate for a full boundary layer simulation; however it provides correct values for the pump performance and allows the analysis in details of the main phenomena involved. The resulting mesh for the entire computational domain under study is shown in figure 3 and some details of the structured mesh are shown in figures 4a and 4b.

Figure 2: Computed flow regime.

Figure 3: Complete grid structure.

Michalis D. Mentzos, Andronicos E. Filios A.E., Dionisios P. Margaris and Dimitrios G. Papanikas

(a)

(b) Figure 4: Mesh views.

3.2 Boundary Conditions The modeled boundary conditions are the ones considered with more physical meaning for turbomachinery flow simulations. At the inlet zone, the axial velocity is a constant based on the through flow for the pump. The absolute tangential velocity at the inlet is zero, which implies in the rotating frame the relative velocity is –rω, and the radial velocity is zero. The involved parameters regarding the turbulence intensity and the hydraulic diameter in the lack of realistic turbulent inflow conditions in industrial applications are estimated with values of 5% and D1/2 respectively. The only specification made at the outlet is that the static pressure in the absolute frame is uniform and is set to zero. This absolute condition is converted into the appropriate relative pressure in the rotating frame. Periodic boundaries are used upstream and downstream of the blade leading and trailing edges, respectively. For the rotating solid surfaces all of the relative velocity components are set to zero, imposing the non-slip condition. 3.3 Numerical solution control The code was run in a 3GHz Pentium IV PC. The number of iterations has been adjusted to reduce the scaled residual below the value of 10-5 which is the criteria. For each run, the observation of the integrated quantities of total pressure, at suction as well as at discharge surface was appointed for the convergence of the solution. In many cases this drives the residuals in lower values than the initially set value. Depending on the case, the convergence was achieved at difference iterations, as the result at a specific mass-flow was used to initialize the computations at another mass-flow. Aiming to smooth convergence, various runs were attempted by varying the under-relaxations factors. In that way a direct control regarding the update of computed variables through iterations was achieved. Initializing with low values for the first iterations steps and observing the progress of the residuals, their values were modified for accelerating the convergence. Since the problem involves both stationary and moving zones, the multiple reference frame model was selected. It is a steady-state approximation in which individual cell zones move at different rotational speeds. As the rotation of the reference frame and the rotation defined via boundary conditions can lead to large complex forces in the flow, calculations may be less stable as the speed of rotation and hence the magnitude of these forces increases. To control this undesirable effect, each run starts with a low rotational speed and then slowly is increased the rotation up to desired level. 4 RESULTS

Indicative results regarding the flowfield simulation of the selected centrifugal impeller, are shown in figures 5, 6 and 7. Pressure and velocity distributions as well as path lines are presented for a volumetric flow rate of Q=65m3/h at n=1450 rpm, which is 60% more than the expected nominal value. The contours of static pressure on the wall of the impeller are shown in figure 5a. It can be seen clearly that the static pressure increases gradually in the stream wise direction and normally the higher pressure is appeared on the pressure surface than on the suction surface on each blade. In figure 5b, the contours of total pressure are shown. It is observed the same quantitatively distribution as in the case of the static pressure. The difference regarding the magnitude of the shown pressures is focused on the term of kinetic energy which is added on the static pressure. The vectors of the absolute velocity close to the impeller walls and to middle of its width are shown in figure 6. Here is clearly visible that the absolute velocity increases as the fluid moves towards to the discharge of the impeller. It is expected since the tangential speed which is the main component varies linearly with the radius. Finally, the path lines colored in agreement with the static pressure variation are shown in figure 7.

Michalis D. Mentzos, Andronicos E. Filios A.E., Dionisios P. Margaris and Dimitrios G. Papanikas

(a)

(b)

Figure 5: Contours of total (a) and static (b) pressure (atm). 2.05e+01

2.05e+01

1.85e+01

1.85e+01

1.64e+01

1.64e+01 1.44e+01

1.44e+01

1.24e+01

1.24e+01

1.03e+01

1.03e+01

8.27e+00

8.27e+00

6.22e+00

6.22e+00

4.18e+00

4.18e+00

2.14e+00

Z

2.14e+00

Z Y

9.38e-02

Y

X

9.38e-02

X

Figure 6: Velocity vectors coloured by velocity magnitude (m/s). 35 1450rpm, 1st order scheme 1450rpm, 2nd order scheme 1950rpm, 1st order scheme 1950rpm, 2nd order scheme

30

25

H (m)

20

15

10

5

0 0

10

20

30

40

50

60

70

80

90

3

Q (m /h)

Figure 7: Pathlines coloured by total pressure (atm).

Figure 8: Predicted H-Q curve at two rotational speeds and for the examined discretization schemes.

Michalis D. Mentzos, Andronicos E. Filios A.E., Dionisios P. Margaris and Dimitrios G. Papanikas

The performance curves of the examined centrifugal impeller for the two studied rotational speeds are shown in figure 8. For each rotational speed, both first order and second order discretization techniques have been utilized. The deviations comparing the applied discretization techniques raises maximum at one percent. In the lack of experimental results the calculations are not validated. However the computation shows qualitatively the trend of the H-Q curves and their relationship applying the similarity rules. 5 CONCLUSIONS

The flow and pressure field through a centrifugal pump impeller results is analyzed solving the fully 3d incompressible Navier-Stokes equations, including the centrifugal force source. Turbulence is simulated with the standard RNG k-ε model. Although grid size is not adequate to investigate local boundary layer variables, global ones are well captured. The present approach is useful for basically understanding the flow states in various operating points. Further work will focus on validating the CFD results with the planned experiments in the near future hopefully at the beginning of autumn. ACKNOWLEDGEMENTS

The authors express their acknowledgements to professor D. Papantonis and lecturer J. Anagnostopoulos in the Hydraulics Machines Laboratory of the NTUA for providing detailed drawings of the examined centrifugal impeller as well as for accessing the experimental facility where the tests of the model pump are accomplished. REFERENCES

[1] Eckardt, D. (1976), “Detailed flow investigations within a high-speed centrifugal compressor impeller”, ASME J. Fluids Eng. 98, pp.390-402. [2] Johnson, M. W. and Moore, J., (1980), “The development of wake flow in a centrifugal impeller”, ASME J. Eng. for Power 102, pp.382-390. [3] Kjork, A. and Lofdahl, L., (1989), “Hot-wire measurement inside a centrifugal impeller”, ASME J. Fluids Eng. 111, pp.363-368. [4] Denton, J.D., (1992), “The calculation of three-dimensional viscous flow through multistage turbomachinery”. ASME J. Turbomachinery 114, pp.18-26. [5] Dawes, M.W., (1992), “Toward improved through flow capability: the use of three-dimensional viscous flow solvers in a multistage environment”, ASME J. Turbomachinery 114, pp.8-17. [6] Casey, M.V., Dalbert, P. and Roth, P., (1992), “The use of 3D viscous flow calculations in the design and analysis of centrifugal compressors”, ASME J. Turbomachinery 114, pp.27-37. [7] Bansod, P. and Rhie, C.M., (1990), “Computation of flow through a centrifugal impeller with tip leakage”, AIAA Paper No 90-2021. [8] Krain, H. and Hoffman, W., (1989), “Verification of an impeller by laser measurement and 3D viscous flow calculations”, ASME Paper 89-GT-150 [9] Farge, T.Z. and Johnson, M.W., (1990), “The effect of backswept blading on the flow in a centrifugal compressor impeller”, ASME Paper 90-GT-231. [10] Zhang, M.J., Gu, C.G. and Miao, Y.M., (1994), “Numerical study of the internal flow field of a centrifugal impeller”, ASME Paper 94-GT-357. [11] Lakshminarayana, B., (1991), “An assessment of computational fluid dynamic techniques in the analysis and design of turbomachinery”, ASME J. Fluids Eng. 113, pp.315-352. [12] Rodi, W., Majumdar, S. and Schonung, B., (1989), “Finite volume methods for two-dimensional incompressible flows with complex boundaries”, Computer Meth. in Applied Mechanics and Engineering 75, pp.369-392. [13] Thakur, S., Wright, J., Shyy, W., Liu, J., Ouyang, H. and Vu, T., (1996), “Development of pressure-based composite multigrid methods for complex fluid flows”. Prog. Aerospace Science, Vol. 32, pp.313-375. [14] Sun, J., and Tsukamoto, H. (2001), “Off-design performance prediction for diffuser pumps. Journal of Power and Energy”, Proceedings of I. Mech. E., Part A, Vol. 215, pp. 191–201. [15] Gonzalez, J., Fernandez, J., Blanco, E., and Santolaria C. (2002), “Numerical simulation of the dynamic effects due to impeller-volute interaction in a centrifugal pump”, ASME Journal of Fluids Engineering, Vol. 124, pp. 348-355. [16] Mentzos, M.D., Filios, A.E., Margaris, D.P. and Papanikas, D.G., (2004), “A numerical simulation of the impeller-volute interaction in a centrifugal pump”, 1st International Conference “From Scientific Computing to Computational Engineering”, Athens, 8-10 September, pp.1-8. [17] Launder, B.E., and Spalding, D.B. (1974), “The numerical computation of turbulent flows”, Complete Methods of Applied Mechanical Engineering, Vol. 3, pp. 269-289.