Numerical investigation of axial force of a centrifugal pump in the reverse mode

This paper numerically investigates the hydraulic performance and axial force of a low specific speed centrifugal pump in the turbine mode using the k - ω SST model. Numerical results are validated against the experimental data in the literature and a good agreement is observed. Results show that the rotation factor increases by flowrate, since the head, which is directly proportional to the transferred circumferential momentum, increases by the flowrate in the turbine mode. Moreover, values of rotation factors on both sides varies in radial direction due to geometrical parameters and, they reduce where larger axial gaps are present. Additionally, the pressure distributions on the hub and shroud sides are correlated with rotation factors on these areas. Interestingly, the axial position of the impeller plays an important role in the axial force of the PAT. By applying impeller movement of 5 mm from its original position, the point of zero axial force has shifted to higher flowrates and the maximum axial force is reduced 50%. Finally, with 10 mm movement of the impeller, the axial force exhibits distinct behavior so that its direction is preserved and, the returned magnitude is almost constant in the working range by 60% decrease in the maximum axial force.


Introduction
Centrifugal pumps can be used for harvesting the stored energy of water streams and supplying to the electrical generators.To achieve this goal, centrifugal pumps are run in the reverse direction and act like turbines.Being inexpensive, available, and easy-to-run, a pump-as-turbine (PAT) can feasibly serve as an alternative to small hydropower plants in the rural and remote areas facing accessibility hardships [1,2,3,4].
Due to the fact that hydraulic characteristic of the pump in the turbine mode are not provided by the manufacturers, choosing the appropriate pump matching the hydropower plant requirements has been an issue [5].The most basic and practical solution to address this problem is predicting the pump performance in the reverse mode using its characteristics in the direct mode.In this method, the PAT characteristics are obtained by the correction factors for discharge, head, and efficiency of the machine in the pumping mode [6,7,8].
In our research group, Derakhshan and Nourbakhsh [9] proposed a correlation for prediction of the whole performance curves of the PAT based on the experimental data of pump with specifics speeds ranging from 14.6 to 55. 6. Comparing this specific speed-based method with the others revealed that it was acceptably valid for prediction of heads and flowrates of the reverse pumps.In a similar research, Barbarelli et al. [10] carried out an experimental investigation on 12 types of pumps in direct and turbine modes to develop polynomial predictive models for pumps running in turbine mode.The proposed method does not present generality, but is straight forward for selecting the appropriate pump for a given hydropower plant.More recently, a novel method for predicting the performance curves and modeling the hill charts of pumps running in turbine mode was proposed by Delgado et al. [11].This method is based on the Hermite polynomial chaos expansion and accurately captures the variable speed operation of a given PAT.
Along with the theoretical methods, Computational Fluid Dynamics (CFD) has been proven to be a valuable tool for estimating the performance of centrifugal pumps and PATs [12,13,14,15,16].To cite an example, Binama et al. [17] carried out an unsteady numerical study on the PAT pressure field under part load using the Ansys-CFX software.They reported that the impeller-volute interaction was the most influential parameter on the PAT pressure pulsation.
In addition to hydraulic characteristics, the prediction of the mechanical forces generated during reverse operation of a pump is vital as they greatly influence the machine life and the power plant sustainability.However, only a few studies concentrated on investigating the mechanical forces in PAT mode in the literature.
In an early work, Gantar [18] experimentally investigated the mechanical forces in an axial PAT.Their findings showed that in the turbine operating mode, the radial force was proportional to the flowrate, and the axial force became about 2-2.5 times greater than the direct mode.For high reliability of the machine, it was recommended to utilize larger bearings for balancing the axial loads.In another investigation, Bario et al. [19], performed a transient numerical study on the radial force of a PAT.They reported that during direct operation, the minimum radial force occurred close to the pump's best efficiency point, while in the reverse mode, the radial force increased with the flowrate.It was observed that in low flowrates, the radial force in turbine operation was smaller than pump mode, while in high flowrates, it became higher.
In our research group, an experimental study on the effect of volute tongue geometry on the performance and radial force of a pump in reverse operation was carried out by Alemi et al. [20].In this study, the tongue with medium stretching was found to reduce the radial force, while more stretching did not return smaller radial force.Also, we showed that the tongue angle could substantially affect both hydraulic performance and radial force.In a follow up research, Alemi et al. [21] attempted to improve the hydraulic and mechanical forces of a PAT by introducing a new volute.Results showed that, for a given impeller, the PAT efficiency in high flowrates is increased by using a volute with a larger cross section, and the radial force noticeably dropped in comparison to the main geometry.
As discussed above, in order to select the most appropriate pump for a given hydropower plant, researchers have developed several methods to predict the hydraulic characteristics of PATs.Meanwhile, mechanical forces affecting the working life and endurance of PATs have not been investigated thoroughly.Only few researches focused on studying the radial forces during the reverse mode, and the axial forces are not yet studied in details.Besides, to the knowledge of the authors, theoretical approaches for describing the behavior of mechanical forces are not included in the literature.Therefore, this research numerically investigates the hydraulic performance and mechanical forces of a centrifugal pump during direct and reverse operation.Meanwhile, the results are to be supported by the physical approaches.
The current paper is organized as follows.In Section 2, the test case and computational details are presented.In Section 3, validations for hydraulic characteristics are discussed.In addition, numerical results for rotation factor, radial distribution of static pressure and axial force are reported.Finally, the main findings of this research are summarized.

Test case
The present numerical study is conducted on the test case presented in [20] wherein experimental data for hydraulic characteristics of a centrifugal pump in turbine mode are reported.This centrifugal pump has six blades with the specific speed of N s = N Q

Governing equations
In order to simulate the fluid flow of the pump, the governing equations of continuity and and 3D steadystate incompressible Reynolds Averaged Navier-Stokes are employed as: where ρ, ν, and −ρu i u j are density, kinematic viscosity and Reynolds stress tensor.To close the above equations, the k − ω SST turbulence model , proposed by Menter [22], is used.This model uses the accuracy of the k − ε model in capturing the bulk flow with the capability of the k − ω model for resolving the near-wall region.
In this model, the Reynolds stress is obtained using Boussinesq eddy viscosity hypothesis (3) as: where δ ij is the Kronecker delta.
One can obtain the turbulent viscosity as: The transport equations for turbulent kinetic energy (k) and specific turbulent dissipation rate (ω) are expressed as: with: where F 1 and F 2 are damping functions.S denotes the invariant of strain tensor, and the constants are as follow: and α 2 =0.44.

Computational grid
The computation domain is composed of six regions namely, impeller, volute, front chamber, rear chamber, inlet and outlet.The fully-developed conditions at impeller inlet and pump outlet are assumed by sufficiently extending the inlet and outlet pipes.In order to generate the computational grid, the commercial software of ANSYS-Meshing was used.Due to the complexity of the geometry of fluid domains, hybrid meshing was implemented.As shown in Figure 2 , in the vicinity of the wall regions, the clustered boundary layer mesh was used, while far from solid walls, the tetrahedral cells were utilized.The grid sensitivity analysis was conducted at the design flow of the pump mode.To this end, four computational grids with the element numbers ranging from 1.40×10 6 to 3.1×10 6 have been generated.

Variations of head coefficient (ψ
)) and hydraulic efficiency of pump (η = ρgHQ/(T ω)) against the number of cells are shown in Figure 3.As it is depicted in Figure 3, the variations of head coefficient and efficiency are negligible when the cell number exceeds 1.7 × 10 6 .Considering the computational costs and the accuracy of the numerical results, the mesh grid with cell number of 1.7 × 10 6 was selected whose details are listed in Table 1.

Numerical simulation
Simulation of the fluid flow is carried out by employing the commercial software of ANSYS-CFX [24].Communication of data between impeller-inlet and impeller-volute domains is performed through the Frozen Rotor (FR) interface while General Connection (GC) is utilized at volute-outlet and volute-front and volute-rear interfaces.The boundary conditions of inlet and outlet are set to total head and constant mass flow respectively.The high resolution advection scheme is used for discretization of the convective terms.This scheme uses second order scheme whenever and wherever possible, and returns to first order when required to achieve a bounded solution.Finally, the criterion for maximum residual of all flow variables is set to 1.0 × 10 −5 .

Results
To facilitate the analysis in the current study, dimensionless flowrate, head and power in the turbine mode are respectively defined as: where φ, ψ and π are flow, head and power coefficients respectively.Ω is the rotational speed and Q, H, T represent the volumetric flowrate, head and hydraulic torque respectively.g and ρ are the gravitational acceleration and water density.b 2 and d 2 stand for the impeller width and diameter at trailing edge.
Accordingly, the hydraulic efficiency of PAT is obtained as: 3.1.PAT hydraulic performance Figure 4 compares the numerical predictions of the pump performance in the reverse mode with the measured data [20] and shows a good agreement between them.As it is shown in Figure 4, the headflow curve of the pump running in the reverse mode follows an increasing trend with the flowrate.It is seen that the head coefficient increases from ψ = 0.16 to ψ = 0.48 by increasing the flowrate from φ/φ n = 60% to φ/φ n = 140%.Figure 4 also displays the power generated by PAT and corresponding efficiency as a function of the flowrate.It can be observed that the power is increasing by the flowrate.As depicted, by increasing the flowrate from φ/φ n = 60%, to φ/φ n = 140%, the power becomes about ten times greater (from π = 0.005 to π = 0.05).Figure 4 also shows the returned efficiency of PAT and compares it with experiments.As it is depicted, the numerical simulations have successfully captured the PAT efficiency with maximum error below 3%.Results show that PAT's efficiency takes the maximum value of approximately 53% and is more sensitive to the flowrate at part load operation.It means that by deviating from the nominal point, the efficiency drops more rapidly when working at low flowrates.Therefore, precautions should be taken so that the PAT operates mostly near its BEP while avoiding part load working.

PAT axial force
Figure 5 illustrates the pump's geometry in the meridional view.In practice, the pressure distributions on the hub and shroud sides of the impeller are not the same, and thus a net force parallel to the rotation axis is generated.Accordingly, investigation of the flow field close to the hub and shroud sides is vital for analyzing the axial force.
The fluid confined within the gaps adjacent to the shroud and hub is not stationary due to the impeller rotation.However, the peripheral velocity in the core flow is almost constant and the confined fluid takes the angular velocity of Γ.For a straightforward analysis, a rotation factor for the core flow is defined as: where k indeed describes the ratio of the core flow angular velocity to the impeller's one.In practice, k is a function of r.However, for convenience, averaged values of k over the shroud and hub sides of the impeller are used.
Under the assumption of the solid body rotation for the fluid in the gaps with the angular velocity of Γ = kω, the pressure distribution on the shroud or hub can be obtained as: Equation (11) implies that as k increases, the pressure distribution in the confined regions declines more rapidly.
Figure 6 displays the rotation factor (k) against dimensionless radial distance (r * = r/r 2 ) on the hub and shroud sides for different flowrates.As it is depicted, k experiences variations near the impeller outer diameter, and thereafter follows a decreasing trend.Next, the values of rotation factors on both sides of hub and shroud increas by flowrate.At the part load operation (φ/φ n = 60%), the rotation factor on shroud is higher than hub, while at off-design operation, hub sides take higher values.The most likely explanations for the behavior of the k are the geometry of confined regions of hubs and shroud as well as the flow regimes in the turbine operation.First, as it is depicted in Figure 5, the shroud region is connected to the volute section by narrow gaps making it more independent from the main stream.Therefore, the changes in the flowrate has less impact on the shroud region than the hub one.Moreover, it can be understood from Figure 6 that in regions where the axial gap between impeller and casing are reduced, the average rotation factor is increased and vice versa.It is due to the fact by increasing the axial clearance, the kinematic energy of impeller due to rotation should be transferred to higher quantity of the fluid, thus the average kinematic energy which is proportional to rotational speed is reduced.This explanation is consistent with findings shown in Figure 6 where k is reduced by decreasing radial distance from r * = 0.7.Moreover, it can be observed that in below r * = 0.6, hub region takes lower values of rotation factor because it has more volume than shroud side in the specified radial distances.In addition, the hydraulic characteristics of PAT (Figure 4) shows that head is increased by the flowrate.Since head (ψ) is directly proportional to the circumferential component of the inlet (c u2 ), it is found that by shifting to higher flowrates, the circumferential velocity of flow in volute, impeller inlet and thus confined regions are increased.Meanwhile, the larger gaps between hub side and volute, the higher rotation factor due to transfer of circumferential momentum.Figures 7 and 8 show the velocity vectors, streamlines and the contour of rotation factor on the shroud and hub sides for the flowrates of φ/φ n = 60% and φ/φ n = 140% respectively.As displayed, in regions near the axis of rotation where the gaps between impeller and casing are increased, the rotation factor is decreased considerably in both side with shroud side taking higher values.Next, the enlarged views of the gaps between hub and casing display that the rotation factor close to the impeller tip is considerably lower than shroud side at φ/φ n = 140% which is related to lower circumferential velocity of the inward mainstream in the turbine operation.However, by shifting to φ/φ n = 140%, the k values are increased considerably due to higher energy and circumferential momentum of upstream flow.
Figure 9 depicts the volume averaged rotation factor over the shroud and hub sides in the reverse operation for different flowrates.It shows that the rotation factor on the hub and shroud sides are increased by flowrate.Although hub side experiences lower values of k in part-load operation, it increases more rapidly than shroud side by the flow and takes higher k values when the flowrate exceeds φ/φ n ∼ = 75%.The rotation factor on shroud commences at k = 0.45 at φ/φ n = 60% and monotonically increases up to k = 0.51 at φ/φ n = 140%.On the hub side, the rotation factor follows more rapid trend in a way that by increasing the flowrate from φ/φ n = 60% to φ/φ n = 140%, the rotation factor rises  from k = 0.44 to k = 0.53.As discussed above, this can be explained by the configuration of impeller so that larger gaps on the hub side cause the rotation factor become more influenced by the flowrate.
Figure 10 gives the radial distribution of dimensionless static pressure (P * = (P − P in )/( 1 2 ρu 2 2 )) over the hub and shroud for the turbine mode.As it is shown, the pressure levels rise with flowrate, which is in line with the previously displayed head-flow curve in Figure 4. Figure 10 shows that by distancing from impeller outer diameter (r * = 1) up to r * = 0.35, the pressure levels are reduced approximately 10% in all flowrates.Next, by increasing the flowrate from φ/φ n = 60% to φ/φ n = 140%, the maximum pressure has become more than two times greater (from P * = 0.29 to P * = 0.66).In accordance with the k values shown in Figure 9, the pressure distributions on the shroud and hub sides become more similar by flowrate up to φ/φ n 80%.However, by moving to higher flowrates, the pressure drop on hub side becomes more rapid, resulting net axial force from shroud to hub side.These findings are supported by the behavior of k on the shroud and hub sides where the intersection of curves occurs close to φ/φ n = ∼ = %75 and moving toward higher flowrate leads to more difference between hub and shroud rotation factors with hub side experiencing higher values.The axial force coefficient of the PAT is defined as: where F ah and F as stand for the axial force on the hub and shroud sides respectively.In addition, A a is the effective area for calculation of the axial force on the shroud or hub side.The coefficient of axial force given by the pressure distributions over impeller sides is also shown in Figure 9.It shows that prior to φ/φ n ∼ = 85% the net axial force is exerted toward shroud side and declines by the flowrate from F * a = 0.5 to zero.Thereafter, the direction of axial force is altered toward hub side and its magnitude is increased to F * a = 1.1 at φ/φ n = 140%.It is worth mentioning that there is a small discrepancy between the flowrates corresponding to zero axial force and the intersection point of rotation factors on the hub and shroud side.This might be due to inevitable numerical errors and the fact that analyzing the averaged rotation factor is an approximation providing insights into the trend of the axial force.Accordingly, the point where k at hub and shroud intersect is not necessarily the same as the point of zero axial force.

Effects of impeller's axial position
As discussed above, different gaps between impeller and confined regions of hub and shroud are considered as the main parameters affecting the behaviors of rotation factor and axial force of the pump in the reverse mode.Therefore, the axial force of the PAT with modified axial positions are investigated.Figure 11 displays the original and modified gaps between impeller and casing after impeller movement of ∆ z = −0.5 mm and ∆ z = −1.0.As it is depicted, minus sign means the impeller has shifted toward hub side which yields larger gaps on the shroud side and smaller ones on the hub side in comparison to the original position.
Figure 12 shows the averaged rotation factor ( k) and normalized axial force (F * a ) of the PAT with axial displacements of ∆ z = −0.5 mm and ∆ z = −1.0mm.As depicted, for ∆ z = −0.5 mm, the hub side takes lower values of rotation factor at low flowrates and increase with higher rate than the shroud side.The k values of hub and shroud sides intersect at flowrate of φ/φ n ∼ = 90%.A comparison between Figure 12 and Figure 9 shows that by moving the impeller ∆ z = −0.5mm, the values of k at hub side approximately remains unchanged.However, the k values on the hub side in increased.Accordingly, intersection point of rotation factors of hub and shroud has moved to higher flowrates.This finding implies that larger gap on the shroud side causes more circumferential momentum transfer from main stream to the shroud region which leads to higher rotation factor in the region of interest.Consistent with the k behavior at ∆ z = −0.5 mm, it can be seen that the point of balanced axial force has shifted to higher flowrates (from φ/φ n ∼ = 75% to φ/φ n ∼ = 108%).In addition, the maximum magnitude of the axial force is reduced by two at high flowrates (from F * a = 1.1 to F * a = 0.5) which is correlated to reduced difference between k values of hub and shroud sides at high flowrates when impeller movement of ∆ z = −0.5 mm is applied.
Figure 12 shows that by further movement of the impeller, ∆ z = −1.0mm, the k at shroud side increase for all working flowrates, while hub side experiences reduction in k.This means that intersection of k values does not occur and therefore, zero axial force is not expected for this impeller configuration.In addition, the difference between k values of hub and shroud regions are almost constant.Consistently, it can be observed from Figure 12 that the axial force for impeller with axial movement of ∆ z = −1.0mm has no direction change, and its magnitude is almost constant over the working range with minor decrease at low flowrates.

Conclusion
This paper presents a numerical study on the behavior of axial force of a centrifugal pump in the turbine mode (PAT).Verification of the numerical simulation is performed by comparing the PAT's numerical results given for head, power and efficiency against the experiments in the literature.Findings of this study reveals that geometrical parameters such as shape of hub and shroud confined regions as well as axial position of the impeller are influential on the resultant axial force.It is found that larger axial gaps between impeller and casing yields lower rotation factor and vice versa.Moreover, the balanced axial force is observed close to the intersection point of rotation factors on hub and shroud sides.Results showed that at the original position of the impeller, zero axial force occurs at φ/φ n ∼ = 80%, and by impeller movement of 5 mm toward the hub side, this point is shifted to φ/φ n ∼ = 108% and, the maximum axial force is reduced by factor of two.Finally, it is interestingly observed that by further movement of the impeller, 10 mm toward the hub side, the axial force experiences no direction change with an approximately constant magnitude.

4 = 10 . 3 .
Its design flowrate and head are of Q n = 8.0 m 3 /h and H n = 13.0 m respectively.The pump has a closed impeller with an outer diameter of 209 mm and operates at the rotational speed of N = 1450 rpm with the Reynolds number of Re = ωr 2 2 /ν = 1.6 × 10 6 .

Figure 5 .
Figure 5. Meridional section of the pump components.The pump configuration imposes larger gaps between the volute and the impeller on the left side (hub) than the right side (shroud).

Figure 6 .
Figure 6.Variation of rotation factor along radial distance at different flowrates. k

Table 1 .
Details of the selected mesh grid.
Interface of Volute with Impeller, Front and Rear Domains Figure 2. Mesh grid.