Analysis of numerical instability of centrifugal pumps operational characteristics

In recent years, computational fluid dynamics has been increasingly used in the development of various types of rotating machines. In the case of water turbines and pumps, a large number of researches have been published recently, related to the use of different numerical methods for prediction of efficiency, cavitation characteristics and different dynamic phenomena. For basic analysis an accurate result near the optimal operating regime is sufficient, but for detailed analysis the numerical analyses in the wider field of operating conditions are necessary. Thus, we encounter some otherwise known physical phenomena that occur in the part load and full load regime. In these areas, due to the nature of the flow, various non-stationary phenomena occur that are dominant for such operating regimes. In this paper, the problems that arises if we consider only stationary results and some recommendations to avoid later problems in the operation of centrifugal pumps are presented.


INTRODUCTION
The pump for a specific application must be selected so that it operates near the best efficiency point (BEP) for as long as possible. The recommended longterm operation of a centrifugal pump is somewhere between 0.6 QBEP and 1.2 QBEP. For the higher specific speed this range is even smaller. This is not always possible and therefore we can get into a situation where the flow conditions are not ideal and different types of vortices and flow separation [1] can form [2]. Such nonstationary phenomena affect the stability of the pump operation and consequently the service life [3] and increased life cycle costs.
With more and more powerful computers and better software, the desire for quality results has increased over time. Comparison of numerical and experimental results showed very good agreement in many cases, but in some cases the results were not acceptable.
In industry where the required time for the development of a new machine is rather limited, development engineers often use the numerical analysis, where they do not consider the time variation of individual characteristics. Due to this fact, the possibility of incorrect simulation is very high. Because it is not always possible to verify numerical results with measurement results, major development errors can be made. These may adversely affect the final result of the development or may extend the development time considerably because more time is needed to devote to experimental work.
The research considers a machine that is very difficult to develop [4] [5]. The reversible pump turbine in pumping mode [6] is analyzed for a large number of different flows and heights. The analysis was also made for different geometric positions of the guide vanes, so that the relative openings of the guide vanes and also the position of the guide vanes relative to the stator vanes in the circumferential direction were varied. It is particularly interesting how the convergence of calculations for each operating point takes place and how large the fluctuations of individual calculated quantities are. In the pumping regime, the flow (Q) was defined in all calculations and the head (H) obtained as a result of analysis. From these calculated results, pump efficiency values can then be determined for all convergence iterations for all specified operating points.
It is important to monitor the convergence of all results and obtain values for all desired quantities in all iterations of the numerical procedure. In some cases, good convergence is observed, where an individual parameter values converges to a solution. In some cases, various size of oscillations is obtained, but are stable and do not increase or decrease. In exceptional situations, the oscillations are constantly increasing and because of too many calculation errors a situation where the numerical calculation diverge may arise.

NUMERICAL METHOD
In this paper, basic relations for the incompressible fluid motion the Reynolds-Averaged Navier-Stokes system of equations; continuity equation and momentum equation are used.
(2) Usual logarithmic relation for the near wall velocity is presented by the equation: , where: , .
(5) The research was based on a multitude of results from different types of centrifugal pumps. For one example, are presented a comparison between stationary and non-stationary calculations and the results of measurements on the model (Fig 1). The comparison shows the trends in the accuracy of individual calculations for the entire operating area. The matching of the results is not the same for all operating points. At higher flows, the differences between the numerical results and the experimental ones are larger. This may be related to the influence of surface roughness [7] [8], which in this case is not considered in numerical calculations. Research has focused on a reversible pump turbine [9]. For the calculation, it was taken a reduced scale model of the reversible pump turbine and operating conditions, defined for the laboratory measurements. The geometry contained a draft tube or diffuser, impeller, distributor, stay vanes and spiral casing. In order to reduce the influence of boundary conditions on the calculation results, at the inlet of the draft tube a fictitious extension that allows the boundary conditions to be as realistic as possible is added, given the current conditions in the various operating modes.
In order to avoid the unintended influence of the computational grids [10] on the result, a quality computational grids with parameters that have been many times verified in other previous cases were prepared. All cases have been confirmed by both numerical and experimental results according to international standards for measurements of hydraulic machines.
In the impeller and the distributor, a structured computational grid (Fig. 2, 3, 4 and 5) with previous checked parameters regarding the values of expansion ratio, aspect ratio and grid skewness were prepared. Depends of the use of the turbulent model, attention was also paid to the size of non-dimensional parameter y + . Based on our previous experience, a k-ω SST turbulent model [11] in all calculations, which has a capability to enables automatic switching between the low-Reynolds model and the use of wall functions was used. A blending function ensures a smooth transition between the two models. An unstructured computational grid in the draft tube (Fig. 6), in stay vanes and in the spiral casing were used, but also followed all the recommendations regarding the quality parameters of the computational grids from a large number of previous numerical analyses. In previous studies, many numerical analyzes were performed, where a comparison was made between the measurement results and the numerical results. For individual types of turbines and pumps and for different operating conditions, it was determined which numerical parameters allow the most accurate results. The boundary conditions in the pumping regime were prepared by defining the flow rate at the outlet of the spiral casing, and at the inlet of the draft tube or the extension of the draft tube, the average value of static pressure was given.
The entire computational grid consists of around 11.000.000 elements:  spiral casing and stay vanes 3.000.000 elements,  guide vanes -3.000.000 elements,  runner -4.000.000 elements and  draft tube -1.000.000 elements.
The computational grids, initial conditions and boundary conditions were the same for all cases of numerical simulations which were compared with each other, so that it was prevented the influence of additional impact on the numerical results.
In generating computational grids, all blades of stator, distributor and runner are considered. This is because, under most operating conditions, the flow conditions in the pumping regime is unsteady and the axisymmetric flow cannot be considered. In non-stationary calculations, the individual time step was equal to approximately 2 degrees of runner rotation, which was shown to be an appropriate size from previous studies. To achieve a stable result in nonstationary analyses, at least three complete runner turns must be analyzed.
Thus, the calculation times for an individual operating point in non-stationary calculations were about 50 hours. These presented average values can differ by ± 20% depending on the position of the operating point in the characteristics. Due to the large number of operating points, two identical workstations were used in parallel to make the research run faster.
When analyzing the flow conditions in the impeller left and right from the BEP point, it can be seen some irregularities in the display of stream lines, but no significant deviations of the flow are observed (Fig. 7). Locally, back flows, flow separation, and vortices can be observed. These are flow phenomena that require suitable turbulent models that are able to describe all the details of the presented phenomena. Above all, these phenomena are such that they change over time and it is therefore necessary to use nonstationary models.
In various operating regimes, certain irregularities in the flow conditions are also found in other parts of the pump. Above all, the input part to the impeller is important, where at flows smaller than optimal the back flow or the input vortex can be noticed.

COMPUTATIONAL RESULTS
Convergence is shown for different geometry configurations and different operating conditions. It can be seen from the results in Fig. 10 that for the BEP, good convergence was actually obtained, so that all the characteristics converge to constant value. Moving in the H_Q characteristic to the left (Fig. 9), some oscillations of the calculated results for different iterations can be observed, although all the computational parameters are the same. In the Fig. 9 the amplitude of efficiency oscillations is around 4%. To achieve the accurate efficiency value, it is necessary to be very careful when to stop the calculation. In the case of flow rate higher than at optimal operating regimes ( Fig. 11 and 12), the magnitude of the oscillations of the calculated parameters is increased to more than 30%. Analyses were performed for several angles of guide vanes opening, and in the paper, the results for angles 19 o and 22 o (Fig. 13, 14, 15 and 16) are presented. In both cases, similar convergence behavior with respect to different operating regimes can be found. The results show that in certain operating points there is a large difference between the results obtained by stationary calculations and non-stationary calculations. Even if stationary results are averaged and compared with non-stationary results (Fig. 17) and (Fig 18), the difference is still almost 3%. In the individual calculations shown as the convergence process, the number of iterations in the graphs differs from case to case, but it is characteristic of all cases that only a smaller number of iterations is shown than was actually calculated for each operating point. Above all, it depends on which operating point was considered and how far from the optimal operating mode it was. They are shown in the main sections, which represent a stable area, even if the amplitudes of the oscillations are very different.

CONCLUSIONS
Precise numerical calculations always require highquality computational grids, which means that the computational area must be divided into a large number of elements. In most cases, this results in long computational times, unless you are using a very powerful supercomputer with a large number of processors. The latter entails a high financial cost, which is not always feasible in the case of industrial research.
All this is already true for stationary calculations, but if the analyses with non-stationary models want to be performed, the computational times are increased enormous so that they are no longer acceptable for industrial research. In industry, more than 90% of Iterations Steady Unsteady Steady_ave numerical analysis is performed using stationary models. In such calculations, however, some problems of oscillations in the convergence of the individual calculated characteristics can be encountered.
In the paper, a problem that is particularly pronounced in numerical analysis of flow in centrifugal pumps, since in this case the flow is decelerated and unsteady phenomena are more likely than in cases where the flow is accelerated, is pointed out. If only the residual convergence is monitored in the calculation, the important information regarding true convergence can be overlooked. Therefore, it is imperative that in the convergence monitoring all the characteristics which are necessary to be obtained by numerical simulation are followed. This allows to obtain information about a potential unsteady phenomenon, which can then be accurately analyzed by an unsteady calculation.
In stationary calculations, the values of oscillations cannot be averaged to obtain the result. It is necessary to obtain results for a large number of time steps and calculate time-averaged values.
The study found that the fluctuations in the results in some operating modes are so large that the result is highly dependent on the moment we stop computing. The difference between the minimum and maximum value can be more than 30%. Therefore, it is good to monitor the convergence of individual characteristics which are important and make an average at appropriate intervals of the number of iterations.
Significantly better results can be obtained from the average values, but the difference from non-stationary calculations is still around 3%. For the first approximation of the result in the development of a new machine, averaging in stationary calculations is approximately acceptable, especially when relative results are sufficient.
For the final result, where the accurate results of the absolute values of certain characteristics are necessary, a non-stationary calculation is recommended, which in most cases means a very long computational time.
Finally, it can be suggested that monitoring the convergence of individual characteristics can be very useful for all operating points, because even in the case when the classical convergence of residues is not the best, certain characteristics can converge even with a small number of iterations. This is due to some local instabilities in the calculation, which, however, do not affect the integral result.
In such a case, the result can be obtained faster and shorten the calculation times, which also means a shorter development time of the new energy machine. This is an important factor of competitiveness in the industry today.