Evaluation of Various Turbulence Models and Large Eddy Simulation for Stall Prediction in a Centrifugal Pump

: Rotational stall is an unstable ﬂ ow phenomenon that reduces the performance of centrifugal pumps, usually occurring under partial load conditions. It causes instability in the ﬂ ow resulting in intense vibrations and noise under certain ﬂ ow conditions. In this study, the one − equation Wray–Agarwal (WA) turbulence model, which was recently developed, is employed to numerically simulate the internal ﬂ ow ﬁ eld of a centrifugal pump under the deep stall condition. The aim of this study is to examine the prediction accuracy for stall by using the WA turbulence model. The method based on computational ﬂ uid dynamics (CFD) has been widely applied for investigation of complex flow patterns in pumps by solving Reynolds − averaged Navier − Stokes (RANS) equations. Particle image velocimetry (PIV) experimental results were compared with simulations predicted using the WA, renormalization group (RNG) k − ε , shear stress transport (SST) k − ω , and realizable k − ε turbulence models and large eddy simulations (LES). The comparisons indicated that the WA turbulence model can accurately predict the flow separation and has a good agreement with the PIV data. The WA model adds a cross − diffusion term and a blending function to the eddy viscosity R equation, so that this model could be expressed as a one − equation k −ω model or one − equation k − ε model as needed by using the switching function. The results show the strong potential of the WA model for accurately computing the stall in rotating fluid machinery. The outcomes of the study are useful in development and optimization of fluid machinery with a low calculation cost and a high accuracy.


Introduction
Centrifugal pumps are among the most commonly used and important general mechanical equipment.When centrifugal pumps are operated at low flow conditions with deviation from the design point, stall is likely to occur due to changes in the inlet angle of attack.Firstly, the flow separation occurs and develops around the suction surface of the impeller blades.Then, the separation vortex blocks the channel and forms a stall cell under deep stall conditions [1,2].The stall phenomenon in centrifugal pumps can be divided into rotating stall and alternative stall.Rotating stall often occurs with an odd number of blades, where the stall cell propagates in the opposite direction to the impellers' rotation [3].The alternative stall often occurs with an even number of blades.The alternative stall refers to the alternating distribution of stall cells in the impeller passages.The circumferential propagation velocity of the stall cell is zero, which is periodically symmetrical [4,5].When stalling occurs in a centrifugal pump, the stall cells destroy the flow field structure with resulting humps and positive slopes in the pump performance curve.It also leads to additional dynamic load on the blades.Under some flow conditions it may cause vibration, fatigue and even fracture of the blades [6][7][8].In addition to reducing the efficiency of the pump, the stall phenomenon also has a serious negative influence on the safety and stability of the pump system [9,10].At present, the research on centrifugal pump stall is still very limited, and therefore in−depth research considering all aspects is needed to explore the mechanism of stall and use this understanding to improve the efficiency and reliability of centrifugal pumps.
Computational fluid dynamics (CFD) is one of the most widely used methods in the design and optimization of fluid machinery [11,12].The advantage of CFD is that it can visualize internal flow fields such as velocity and pressure.It is also useful to optimize the geometry parameters to enhance the hydraulic performance of fluid machinery [13].Turbulence models play an important role in achieving accurate and efficient simulation results using RANS equations.The direct numerical simulation (DNS) method employs the transient Navier−Stokes equations to solve the turbulent flow and in principle can provide very accurate results [14].However, DNS requires very small spatial grids and step sizes to resolve the detailed spatial structure and temporal properties of turbulent flow, which can be very computationally intensive and may not be feasible on the currently available computing hardware [15] for the flow conditions in fluid machinery.The flow in fluid machinery is turbulent with high Reynolds numbers, where large scale vortex structures play the leading role and small−scale vortex structures play the role of transmission and dissipation of energy.Based on this idea, Smagorinsky proposed the large eddy simulation (LES) method in 1963.The fundamental principle of LES is to filter the small−scale eddy structures through a filtering function and then directly calculate the large−scale eddy structures and their effects on the turbulent flow solution; the other small−scale eddies are modeled using a sub-grid-scale Smagorinsky model [16].LES is less computationally intensive compared to DNS but still requires a higher near−wall mesh density than the RANS solutions and therefore more computing resources.Until now, RANS is the most commonly used method for the numerical simulation of fluid machinery [17,18].
RANS equations do not solve the instantaneous Navie−r−Stokes equations directly, but average the transient fluctuations and employ a turbulence model to model the turbulent stresses and then solve the simpler time−averaged equations [18].In the past five decades, various RANS turbulence models have been developed.The difference among these turbulence models lies in the mathematical modeling of Reynolds stresses employing different assumptions [19].Depending on the number of partial differential equations used to calculate the eddy viscosity, turbulence models are usually classified as zero−, one− and two−equation models [20,21].At present, the two−equation turbulence models are the most widely used for prediction of flow field characteristics of pumps under design flow conditions [22].However, when the flow deviates from the design flow condition, the vortex viscosity assumption will distort the calculated average flow due to the separated flow and the influence of vortex flow.The simulation results of different models under off−design conditions are different [23].Furthermore, the prediction of complex 3D flow needs a large number of grids to solve the RANS equations, which requires larger computational resources.Therefore, it is important to use fewer transport equations to develop more accurate and efficient turbulence models [24].One−equation turbulence models are used to calculate eddy viscosity directly via solving the transport equation for eddy viscosity rather than having to calculate two or more turbulence variables.One−equation turbulence models are simpler and require less computational time compared with two−equation turbulence models.However, various existing one−equation turbulence models have low accuracy in predicting non−equilibrium flows [25].Therefore, it is important to improve the prediction accuracy of one−equation turbulence models for a better understanding of turbulence.
The recently developed Wray−Agarwal (WA) turbulence model is a one−equation model based on k−ω closure [26] and shows great potential in predicting complex three−dimensional (3D) flows [27][28][29].The WA model showed higher accuracy than the Spalart−Allmaras (SA) model for many flows, and can accurately simulate non−equilibrium and low Reynolds number effects in the near−wall region.Han et al. [30] used the WA model and four other commonly used one− or two−equation turbulence models to predict the internal flow of a three−dimensional U−bend.The results showed that the velocity iso−surface calculated using the WA model is largely consistent with the test results.Ouyang et al. [31] employed the WA model to calculate the multiphase flow in a rotating packed bed.Their results showed that the WA model is superior to the other two−equation models in predicting free shear flow.Ji et al. [24] used the WA turbulence model for a centrifugal pump with semi−spiral suction chambers, the results showed that this model could calculate the energy performance of centrifugal pumps effectively in the entire operating range, and had higher accuracy than other models.In over 70 applications to date, the WA model has shown its ability to predict complex 3D flows quite accurately and, in many instances, exhibited superior performance compared to SA and other two−equation models based on k−ε closure.
The operation of centrifugal pumps deviates from the design conditions, resulting in a very complex internal flow field, which has high requirements on the efficacy of a turbulence model.Pederson [32] used LES to predict the flow field inside an impeller at 0.25 times the design operating conditions.The existence of the alternating stall phenomenon and associated vortices was successfully verified by comparison with test results.Zhang et al. [33] used k−ε and SST k−ω turbulence models to simulate the flow field of an impeller under off−design conditions.The two turbulence models can get accurate results overall, but there were some deviations locally from the experimental results.Therefore, improving the simulation accuracy of turbulence models for the application of centrifugal pumps is a current research hotspot.
In this study, the WA turbulence model and the realizable k−ε, the RNG k−ε, the SST k−ω and LES are employed to simulate the impeller passage of a centrifugal pump under strong part−loading conditions firstly.Then, the numerical simulation results are compared with the PIV test data [4].The purpose of this paper is to compare and validate the accuracy of the WA turbulence model in predicting stall phenomena compared to other turbulence models.

Wray−Agarwal (WA) Turbulence Model
The innovation of the WA model is that it includes a cross−diffusion term and a blending function in the eddy viscosity R (=k/ω) equation, such that the two destruction terms can be switched smoothly [26,34].The transport equation of the variable R is given by: The equation for the average strain S is Equation ( 3) is the turbulent eddy viscosity equation: The blocking effect of the wall can be explained by the damping function fµ: By employing the switching function f1, the WA model can be expressed as a one−equation k−ω model in the near−wall region and as one−equation k−ε model away from the wall, as listed in Table 1.For higher stability, the value 0.9 is set as the upper limit in the switching function f1:  = min tanh arg , 0.9 (5)

Large Eddy Simulation (LES) Method
The large−scale fluctuations that are directly calculated are called solvable or resolvable turbulence scales, and the small−scale fluctuations that are not directly calculated are called the non−resolvable−scale or subgrid−scale turbulence [35,36].
The filtered equations are obtained by filtering the Navie−r-Stokes equation: Here,  is the component of the filtered velocity in the , ,  direction, and  is the filtered pressure. is the subgrid−scale (SGS) stress, which reflects the effect of small−scale motion on large−scale motion and is defined as  =   −   .
The SGS stress tensor due to filtering is unknown and needs to be modeled.The specific expressions used in the modeling are given below: Germano et al. [35] proposed a dynamic Smagorinsky model (DSM), which was followed by Lilly [36] improving the model.The basic idea of a dynamic model is to introduce the local structure information of turbulence into SGS stress through two filtering steps, and then adjust the modal coefficients in the calculation process.The stress tensor can be expressed as The tilde represents the filtering operation with scale ∆ ;  is the resolved velocity component

Computational Model
The model selected is a multistage centrifugal pump commonly used in the industry, which was tested in the Fluid Mechanics Laboratory of the Technical University of Denmark [4,32].PIV technology was used to study the flow field inside the impeller under the deep stall condition.After the release of the PIV experimental data, this model has been investigated by many scholars who have verified the reliability of the experimental data.The impeller structure is shown in Figure 1.The main geometrical data and design parameters of the impeller are shown in Table 2 [37].

Meshing and Boundary Conditions
Figure 2a shows that in this study the TurboGrid 2019 R3 software was used to generate a high−quality hexahedral mesh and that the mesh around the blades was refined.Since the impeller domain has a circumferential symmetry and six blades, a single channel is meshed and then replicated six times in the circumferential direction.Mesh nodes connected to replication domains merge smoothly at the domain boundaries.The value of y + in the whole impeller domain is less than 1, as shown in Figure 2b.

Mesh Independence Verification
The choice of mesh size should be such that the mesh density has no significant effect on the numerical simulation accuracy.Four mesh sizes (2.4 × 10 6 , 4.3 × 10 6 , 5.4 × 10 6 , 7.4 × 10 6 ) were chosen to calculate the sensitivity of the grid number to the numerical simulation results.A value of 0.25Qd shows a strong part−load and unstable operating condition, the head and efficiency of the pump may fluctuate significantly due to the influence of internal vortices.Therefore, the method of obtaining pressure data by setting monitoring points located inside the impeller will make the results more accurate.As shown in Figure 3a, six monitoring points are selected in the middle section of the impeller from the cross−sectional line of the radial position r/R2 = 0.6.The differences in the numerical simulation results of the pressure values at each point for different numbers of meshes were analyzed.Figure 3b shows the pressure values at each point with different mesh numbers.When the grids have small numbers, the pressure values at the same point fluctuate considerably.But when the mesh number reaches 5.4 × 10 6 , the change in the pressure value is less than 4%.It shows that the flow field can be accurately predicted when the mesh number is 5.4 × 10 6 .To ensure accuracy and efficiency, the mesh size 5.4 × 10 6 was chosen for the following numerical calculations.To reach computational stability as soon as possible, the results of solving a steady flow field with the SST k−ω model were used as the initial flow field for the LES model calculation.The calculation results were recorded after 10 rotations of the impeller.The SIM-PLEC algorithm was chosen to couple velocity and pressure.In all simulations, the residuals of the continuity and momentum equations, and the convergence criteria of the centrifugal pump head were set as 10 −5 .

Experimental Model
In the experiments at the Technical University of Denmark [4,32] PIV experiments were performed using CCD cameras, particle imaging velocimeters, pulsed laser emitters, image converters and synchronizers.PIV was performed on the horizontal (x, y) plane inside the impeller, which was located 2.9 mm above the plane hub, and on half of the impeller outlet width (Z/b2 = 0.5).The locations of the two 93 × 94 mm 2 observation windows are shown in Figure 4.The relative velocity of different radial position r/R2 at the impeller's middle height was studied.As shown in Figure 6, the velocity value on the suction side in channel A predicted by the PIV is slightly larger than that on the pressure side.The reasons are twofold: one is that the stall cell at the entrance of channel B pushes a large amount of fluid into channel A; the other is that the internal flow is superimposed by the axial flow and a uniform flow [38].The realizable k−ε, RNG k−ε, SST k−ω and LES models predict the presence of larger stall cells in channel A, with higher relative velocity values on the pressure side compared to the PIV experiment.In addition, the distribution of relative velocity values in channel B predicted by the four models, except for the WA model, are irregular.The relative velocity calculated using the WA model can match the experimental results, especially the simulation of the non−stall channel is in much better agreement with the experiment than the other models are.Here, the velocity distribution obtained with the four turbulence models under a 0.25Qd condition, the radial and tangential velocities at r/R2 = 0.50 and 0.90 for each turbulence model, are analyzed.The positions of selected monitoring lines are as shown in Figure 7.As shown in Figure 8, the velocity value on the suction side of channel A inlet is much higher than that on the pressure side.However, the calculation results of the three models realizable k−ε, RNG k−ε and SST k−ω are similar to each other but are contrary to the experiment.This significant difference is mainly attributed to the effect of inlet pre−rotation [37].Under a small flow rate condition, the backflow generated by the rotation of the impeller combines with the leakage flow resulting in a large number of eddy currents at the inlet, and the influence of the eddy currents is ignored in the calculation process of these three models.The WA model considers the influence of the inlet eddy current, and the velocity distribution is very close to the experimental results.In channel B, all five models and experimentally obtained velocity maps show stall phenomena.The RNG k−ε, SST k−ω, LES and WA model results have good agreement with the PIV results with slightly different peaks.

Internal Flow Field Analysis
The flow velocity at the exit of channel A is more uniform than at the entrance.The stall cell predicted by the realizable k−ε, RNG k−ε and SST k−ω models in channel A do not cover the exit; thus, the five models have achieved good prediction results.At the outlet of channel B, the radial velocities calculated using the SST k−ω, WA and LES models coincide well with PIV test in trend.In general, the LES and WA model results are closest to the PIV data for the velocity prediction in the entire impeller channel.

Vorticity and Vortex Intensity Analysis
Next, the five turbulence models are evaluated by computing the Z−directional vorticity magnitude.Figure 9 is a graph showing the intensity distribution of the Z−direction vorticity component.From Figure 9f [32] it can be observed that small stall cells shed from the impellers' leading edge in channel A are distributed near the suction side of the blade and migrate downstream.In channel B, a small clockwise rotating vortex is shed from the pressure leading edge, and the area with high vorticity occupies the space of the impeller's inlet.The five turbulence models successfully predict the two regions with high vorticity values in the non−stall channel, but the eddy rotation direction of the five models is opposite to that of the test in some ranges.The realizable k−ε model overestimates the two regions with high vorticity values and has the lowest agreement with the experiment.The predictions of the RNG k−ε, SST k−ω and LES models for the high vorticity region in channel B match well with the PIV data, but overestimate the range of stall vortices in channel A. The WA model has good agreement between the predicted vorticity value and the experiment, and the initial point and range of the stall cell are predicted successfully.The distribution of the vortex structure in the stall cell was analyzed using the Q criterion [37].The Q criterion can be defined as: where Ωij is the rotation rate tensor, Sij is the strain rate tensor, and cq is the empirical coefficient (usually 1).When Ωij > Sij, the rotational force on the fluid is greater than the stress force, and it is considered that there is a vortex structure present.The vortex intensity increases with an increase in Qc.
The stall cells are analyzed using the Qc iso−surface method.The value of the iso−surface has great influence on the expression of the stall cells.In order to express the internal vortex structure intuitively and clearly, a stall vortex of Qc = 3000 1/s 2 was analyzed.The vortex structure and streamline diagram inside the impeller are shown in Figure 10.There are many stall cells gathering at the blade outlet, and the main reason for these vortex structures is the wake jet.For the realizable k−ε, RNG k−ε and SST k−ω models, the stall cells in stall channels are significantly more than those in non−stall channels.LES model recognition of vortex structures on the same iso−surface is more abundant, much more than by the other models.The inlet reflux vortex simulated by the WA model is located above the stall channels, while the simulation results of the other models are located above the non−stall channels.

Figure 2 .
Figure 2. Mesh of the impeller and y+ distribution on the blades.(a) Structured mesh in the impeller; (b) y + distribution.

Figure 3 .
Figure 3. Mesh independence of solution analysis.(a) Position of monitoring point; (b) pressure changes with different mesh numbers.The inlet boundary condition was set as total pressure, and the outlet boundary conditions was set as mass flow.The entire computational domain was discretized by the finite volume.Non−slip wall was set as the wall boundary condition.RNG k−ε, realizable k−ε, SST k−ω and LES are already available in ANSYS−Fluent.By employing the User Defined Function (UDF), the WA turbulence model was imported into ANSYS−Fluent.The time step in LES model calculation was set as 3.5 × 10 −6 s, and the maximum iterations in each time step were 20, so as to ensure that the courant number was controlled within 1.To reach computational stability as soon as possible, the results of solving a steady flow field with the SST k−ω model were used as the initial flow field for the LES model calculation.The calculation results were recorded after 10 rotations of the impeller.The SIM-PLEC algorithm was chosen to couple velocity and pressure.In all simulations, the residuals of the continuity and momentum equations, and the convergence criteria of the centrifugal pump head were set as 10 −5 .

Figure 5 Figure 5 .
Figure 5 is the streamline distribution on the impeller's middle section using all four turbulence models and PIV experiments under a 0.25Qd flow condition.Channel A is the non−stall channel and channel B is the stall channel.Under 0.25Qd, the non−stall channel and stall channel in the impeller are distributed alternately.The flow in channel A is smooth and regular, and the streamlines are consistent with the blade profile.There are only a few mild flow separations around the blade surface.In channel B, there are two large stall cells.The stall cell can strongly block the inlet of this channel, resulting in the reduction in the flow capacity of this channel.Around the blade's pressure surface, the stall cells occupy most of the area of the channel, and thereby push the fluid to the suction surface.All models predict the alternating stall phenomenon successfully; however, they all show different flow field details.The realizable k−ε, RNG k−ε and SST k−ω models predict that the suction surface of channel A has a large stall cell, which is significantly different from the experimental results.The realizable k−ε model considerably overestimates the stall cell size at the entrance of channel B. The RNG k−ε and SST k−ω models overpredict the number of stall cells in channel B, but the streamlines in channel B have good agreement with the PIV results.The LES predicts smaller stalled cells on channel A, which is a small difference from the PIV test results.The LES model predicts flow in channel B in good agreement with the test data.The WA model predicts flow patterns very well in channel A, reaching almost the same values as in the PIV experiment.The number of stall cells in channel B and the separation flow point obtained from the WA model coincide well with the experiment, but the range of stall cells is slightly different from the experiment.In general, the LES and WA models have good agreement with the experiments.The WA model could accurately calculate the flow separation in the impeller under an adverse pressure gradient.

Figure 7 .
Figure 7.The positions of selected monitoring lines.

Table 1 .
Model constants of the WA model

Table 2 .
Impeller geometry and design parameters.