Numerical Investigations of Unsteady Flow in a Centrifugal Pump with a Vaned Diffuser

Chalmers Publication Library (CPL) offers the possibility of retrieving research publications produced at Chalmers University of Technology. It covers all types of publications: articles, dissertations, licentiate theses, masters theses, conference papers, reports etc. Since 2006 it is the official tool for Chalmers official publication statistics. To ensure that Chalmers research results are disseminated as widely as possible, an Open Access Policy has been adopted. The CPL service is administrated and maintained by Chalmers Library.


Introduction
In centrifugal pumps, the relative motion between the rotor and stator and the small radial gap between the impeller blades and diffuser vanes result in a highly unsteady flow.This unsteadiness creates high pressure fluctuations, which are in turn responsible for unsteady dynamic forces that create vibrations and can cause damage.A large amount of detailed experimental investigations have therefore been dedicated to understanding the flow in centrifugal pumps.Among those, Dring et al. [1] showed that the two major sources of unsteadiness are potential and blade/wake interactions.In centrifugal turbomachines, the effects of these sources of unsteadiness become comparable [2].On the basis of the studies mentioned above, Ubaldi et al. [3] built a simplified model of a centrifugal pump with a rotatable vaned diffuser to study rotor-stator interaction.They then investigated the upstream effect of the vaned diffuser on the impeller outflow in the radial gap of the model, as well as the flow in the impeller [4][5][6].
The experimental work contributes to an understanding of the flow complexity owing to rotor-stator interaction in the centrifugal pump.However, the knowledge is limited to the number of measurement points.For an extensive and detailed analysis of the flow, many different probes must be positioned in the geometry, although the complete flow field is not monitored.Computational fluid dynamics (CFD) techniques have been shown in the recent decades to be a useful complement to experiments.CFD calculations can provide more extensive results in the whole domain, giving a better overall understanding of the flow in the whole turbomachine.In recent years, improved computational algorithms and hardware development have shown convincing evidence that CFD calculations are reliable tools that can be used to analyze the unsteadiness of the flow [7].However, the methods and software used to make the CFD calculations must be validated by experiments.To achieve this, the European Research Community on Flow Turbulence and Combustion (ERCOFTAC), together with Ubaldi et al. [3], adopted the centrifugal test rig as a test case for joint experimental and  theoretical investigations of rotor flow and rotor-stator interaction.The original test case was presented by Combes [8] at the Turbomachinery Flow Prediction ERCOFTAC Workshop in 1999.Intensive studies were then carried out using proprietary CFD software.2D numerical simulation was initially the only approach permitted by computer hardware limitations.Bert et al. [9] presented a 2D analysis of the ERCOFTAC centrifugal pump.Following the development of hardware, 3D unsteady studies were then done [10,11].The geometry of the ERCOFTAC centrifugal pump is shown in Figure 1.Large meshes and short time steps are often used in 3D unsteady calculations, making the simulations computationally heavy.Simulations of this kind are decomposed for parallel processing.This becomes costly when proprietary software is used, where there is an additional license cost for each process.To offer a viable alternative, the community-driven OpenFOAM Turbomachinery Working Group extends and validates OpenFOAM for turbomachinery applications [12].
OpenFOAM is an open source library written in C++ [13].It is based on the finite volume method and has proven to be as accurate as proprietary codes for many applications [14][15][16].2D numerical simulations of the ERCOFTAC centrifugal pump were previously made using OpenFOAM by Petit et al. [17].The general behavior of the flow was well captured, but the results suggested the use of 3D simulations for better capturing the unsteadiness of the flow.
The present work reports the unsteady flow field of the ERCOFTAC centrifugal pump obtained by 3D steady and unsteady CFD calculation.The steady simulation uses the frozen rotor concept, where the results are a crude estimation of the ensemble-averaged flow for a fixed rotor position.A series of such snapshots gives an estimation of the unsteadiness of the flow in the pump.The unsteady simulation uses a fully resolved sliding grid approach.The unsteady flow is computed using four different turbulent models, the  − , the realizable  − , the RNG  − , and the  −  SST.The results are analyzed and compared in detail with the measurements performed by Ubaldi et al. [3] in the radial gap between the impeller and diffuser.To this day, the present work is the most extensive and accurate comparison between experimental and numerical results of the ERCOFTAC centrifugal pump unsteady flow field.All the available experimental results are compared with the numerical results, to analyse the accuracy of the two main approaches used to simulate rotor-stator interaction.It is furthermore an open test case, that has been shared with the OpenFOAM community, and the results presented in the present work can be easily reproduced.

Test Case and Operating Conditions
The ERCOFTAC centrifugal pump test rig was built by Ubaldi et al. [3] and consists of a 420 mm diameter unshrouded centrifugal impeller and a 644 mm diameter radial vaned diffuser.Details on the geometry and coordinates of the impeller blade and diffuser vane profiles are given in Ubaldi et al. [3].The impeller has 7 untwisted constant thickness   backswept blades, and the diffuser has 12 vanes.There is a 6% vaneless radial gap between the impeller and diffuser.The tip clearance can be varied but is set at a value of 0.4 mm, corresponding to 1% of the blade span.The tip clearance is not included in the present simulations however.
The measuring techniques used were hot wire anemometry and fast-response pressure transducers.The hot wire probe was used to measure the unsteady 3D flow in the vaneless gap at a radial distance of 4 mm from the trailing edge of the blade and 8 mm from the leading edge of the vane.The unsteady static pressure was measured at the front cover facing the unshrouded impeller passages.The experiments were conducted at a constant rotational speed of 2000 rpm, at the nominal operating condition described in Table 1.The model operates in an open air circuit directly discharged into the atmosphere from the radial diffuser.The inlet air temperature was 298 K and the air density was 1.2 kg/m 3 .

Computational Domain and Numerical Setup
The meshes for the impeller and diffuser regions were generated separately; see Figure 2.Both of the mesh regions are block structured and were generated using the ICEM-CFD software, with O-grids around the blades and the diffuser vanes and H-grids in the blade passages.Angular geometrical periodicity is used to mesh the different blades.The mesh has    turbulence Reynolds stress tensor, and the standard log-law wall function is used at the walls.A second-order, linearupwind scheme [18] is used to discretize the convection terms, and the second-order backward scheme is used to discretize the time derivative.The number of iterations in the transient SIMPLE algorithm at each time step is set to 10.This number of iterations is enough to reduce the residuals by three orders of magnitude.The final residuals are below 10 −6 except for the pressure, which has a final residual of 10 −5 .Based on the experimental measurements [3], the inlet velocity is set to a constant purely axial component, and the inlet turbulence intensity is assumed to be 5%, with a viscosity ratio of 10.At the outlet, all variables are given a zero gradient boundary condition, except for the pressure, which is set to 0.   The time step for the unsteady simulations was fixed to 7 −5 s, which gives a maximum Courant number of 6.The calculations are carried out on a cluster, using 8 nodes equipped with 4 cores each.Using these parameters, one impeller rotation in the unsteady simulation takes 48 hours to compute.
Although the impeller tip clearance is a very important aspect of the geometry, it is not included in the computational model.The tip gap is 1% of the blade spade, that is to say, 0.4 mm.To validate the log-law wall function used at the walls, the  + values at the walls must yield a minimum value of 30.This in turn gives a limitation for the cells size at the walls bigger than 1% of the blade spade.It is thus very difficult to include the tip clearance in the computational domain when wall treatment is applied at the walls.
Two different rotor-stator interaction methods are used in this work.The first is a steady-state multiple frame of reference approach, where the impeller flow field is solved in a rotating frame of reference and the diffuser flow in a fixed frame of [17].The local fluxes at the interface are transferred using a General Grid Interface (GGI), implemented in OpenFOAM by Beaudoin and Jasak [19] and validated by 2D  simulations of the ERCOFTAC centrifugal pump [17].It is a method commonly used in proprietary CFD codes, where it is referred to as frozen rotor.This approach resembles a snapshot of the flow, frozen in time, thus failing to predict the unsteadiness of the flow between the diffuser vanes.It also gives an incorrect advection of the impeller wakes in the vaned diffuser.It is a fast preliminary method, however, and gives satisfying results in the radial gap between the impeller and diffuser [20] that can also serve as a good initial condition for unsteady simulations.The second rotor-stator interaction method resolves the flow unsteadiness using a sliding grid approach.This transient method rotates the impeller part of the mesh with respect to the stator part at each time step.In this case, the local fluxes are also transferred using the GGI interface, which is updated every time step.The interaction between the impeller and diffuser is thus fully resolved.The chosen time step corresponds to a rotation of the impeller of about 0.84 ∘ .Previous studies have shown that this angular time step is sufficient to allow the transient solver to catch the unsteadiness of the flow [21].The results from two different rotor-stator interaction approaches are first compared with each other and with the measurement database, using the standard  −  turbulence model to model the eddy-viscosity.A comparison is then made of the results from four different turbulence models.An extensive comparison of the different RANS turbulence models available in OpenFOAM was made by Moradnia [22].The focus in Moradnia's study was particularly on the highand low-Reynolds-number two-equation eddy-viscosity turbulence models and on the Reynolds stress transport turbulence models.The conclusion of the study by Moradnia was that the four high Reynolds turbulence models were better at predicting the flow while preserving numerical stability and efficiency.These are the  −  [23], the realizable  −  [24], the RNG  −  [25], and the  −  SST [26] turbulence models.The  −  SST is used in combination with the standard loglaw wall function at the walls.These models were thus been chosen for comparisons in the present work.

Results and Discussion
The following sections present and discuss the results of the steady and unsteady rotor-stator interaction approaches using the standard  −  turbulence model.The results are compared with the experimental database.We also discuss the main differences between the results from the frozen rotor and sliding grid approaches.Finally, comparisons of the unsteady results using four turbulence models,  − , realizable  − , RNG  − , and  −  SST, are presented.
All the kinematic quantities are normalized with the rotor tip speed,  2 .The circumferential rotor relative coordinate,   , and the axial coordinate, , are made nondimensional by means of the circumferential pitch of the rotor,   = 2/  , and the blade span at the rotor outlet, , respectively.The velocity comparisons are made at midspan in the radial gap, at a radius of 1.02 *  2 .The time, , is normalized by the rotor blade passing period,   .Triangles and squares in the plots describe the positions of the impeller blades and diffuser vanes, respectively.The results are observed in the relative frame, such that the impeller blades (triangles) have a fixed position in time, while the diffuser vanes (square) move from left to right.

Steady-State Results with the 𝑘−𝜀 Turbulence Model. Figures 3 and 4
show the ensemble-averaged experimental and calculated radial and tangential velocities in the radial gap at  = 1.02 *  2 over two impeller pitches between the hub (/ = 0.05) and the shroud (/ = 0.95).The behavior of the flow is roughly captured by the simulation.However, a shift of the predicted wakes can be observed in Figure 4. Figure 5 extracts the same results at mid-span (/ = 0.5).Although the radial velocity is reasonably well predicted, the phase shift can also be observed.The tangential velocity is incorrectly predicted after the rotor blade passage, as the wake remains at the same place, and aligns with the local relative streamlines, instead of being adverted into the diffuser vanes channel and dissolving.
/  represents the position of the impeller blade relative to the diffuser vane.For steady simulation, /  does not vary over time.The rotor part of the mesh must be rotated to get a new position for the impeller blade, and the flow must be computed once again.A sequence of many steady simulations with different angular positions for the impeller blades will then give a sequence of snapshots of the unsteadiness of the flow.using the  −  turbulence model are presented in Figures 6  and 7, respectively.The calculated profiles are in good agreement with the measurements.They show a good prediction of the wake position and of the potential interaction linked to the diffuser vane leading edge, corresponding to low radial velocity peaks.However, some differences are still visible.The calculations overpredict the maximum of the wake for the tangential velocity, and a slight phase shift is still visible at the maximum of the wakes.This can be explained by Figures 8 and 9, where it can be seen that the wakes predicted by the calculation are orthogonal to the shroud, while the measured ones are slightly inclined.This creates a phase shift at mid-span.A comparison of the static pressure is presented in Figure 10.The general evolution of the pressure is well predicted, but a slight phase shift is also present, and the global fluctuation of the mean pressure level is not accurately predicted.It is important to remember that some secondary flows were not included in this calculation.For example, Ubaldi et al. [6] explained that the low-pressure region located in the blade passage near the suction side seemed to be caused by the tip clearance vortex.

Unsteady Results with the
Figure 11(a) shows the pressure signals obtained by two probes located in the radial gap and at the outlet of the diffuser.The frequency contents of the signals are shown in Figure 11(b).The reconstructed signals using selected frequencies are shown in Figure 11(c).The frequency of the rotor is   = 2000/60 = 33.3Hz.The frequency of the impeller blades is   =   *   = 232 Hz and that of the diffuser vanes is   =   *   = 400 Hz.The signal in the radial gap is dominated by the impeller blades' frequency and its multiples,   ; see Figure 11(b).The impeller wakes are high-energy regions and decay very rapidly in the radial gap.At the outlet of the diffuser, a similar signal is shown in Figure 11(c), but a reduction in amplitudes with values lower than those found in the pump is apparent.

Comparison of the Frozen Rotor and Sliding Grid
Approaches. Figure 12 illustrates the main difference between the frozen rotor and sliding grid results.The velocity magnitude behind the rotor blades in the frozen rotor simulation is shown in Figure 12 problem with the frozen rotor concept.The wakes advected in the sliding grid approach do not change direction when they enter the diffuser region but are advected in the diffuser vane passages and cut by the diffuser vanes.

Comparison between Four Turbulence Models.
The four two-equation turbulence models,  −  [23], realizable  −  [24], RNG  −  [25], and  −  SST [26], used with the sliding grid approach are here compared for one position of the impeller blades with respect to the diffuser.The position differs for the pressure and velocity comparisons since Ubaldi et al. presented the velocity measurements for /  = 0.126 and the pressure measurements for /  = 0. Figures 13 and 14 show that all four turbulence models predict the flow distribution very well.The differences between the numerical and experimental results may be explained by the lack of tip clearance secondary flow in the numerical model, as suggested by Ubaldi et al. [3].It is interesting to note that the calculation using the realizable  −  turbulence model very accurately predicts the maximum and minimum peak values of the wakes for the radial velocity, although the wakes are thinner than the ones of the experimental results.However, the tangential velocity in the impeller wakes is largely overpredicted when the realizable  −  turbulence model is used.Figure 15 shows that all four turbulence models predict the pressure drop in the wakes behind the impeller but fail to predict the pressure fluctuations in the diffuser vanes passage.In general, the four different turbulence models seem to overpredict the tangential velocity, as shown in Figure 16.This could be explained by the lack of leakage flow, as well as the use of the wall function approach, which fails to predict an accurate boundary layer flow.Figure 17 shows that all the models predict the large radial velocity at mid-span near the pressure side of the blades.This was not predicted at all in previous studies by Combès et al. [27] and Petit et al. [17].The calculations made with the  −  and the RNG  −  turbulence models predict this feature in the lower part of the channel rather than that centered in the mid-span, as shown in Figure 17.

Conclusion
3D steady and unsteady simulations were carried out with the OpenFOAM CFD tool to study rotor-stator interaction in a centrifugal pump with a vaned diffuser.The four basic turbulence models,  − , realizable  − , RNG  − , and − SST, were used to perform calculations.The results were compared with the experimental investigation presented by Ubaldi et al. [3].Although the secondary leakage flow was not included in the numerical model, the results show qualitatively good agreement with the experimental results.
The steady numerical simulation with the − turbulence model shows reasonably good agreement with the measurements, although it obviously does not predict the important unsteadiness.The numerical results in the radial gap yield good agreement with the experimental ones, but the flow in the diffuser passage is incorrectly predicted.The corresponding unsteady simulation resolves most of the unsteadiness and accurately predicts the behavior of the flow.Analyzing the pressure signals at the impeller and diffuser exits, it was observed that the impeller wakes are advected to the exit of the diffuser.
The results of four different turbulence models were compared with the experimental results.The turbulent models yield slightly different results but predict the flow unsteadiness at the same level of accuracy.
To further improve the numerical results, the leakage flow between the impeller blades and the shroud should be taken into account.Large eddy simulation (LES) or hybrid methods such as detached eddy simulation (DES) are suggested to more accurately resolve the flow unsteadiness.

Figure 3 :
Figure 3: Experimental and calculated radial velocities in the radial gap for the steady-state simulation using the  −  turbulence model.

Figure 4 :
Figure 4: Experimental and calculated tangential velocity fields in the radial gap for the steady-state simulation using the  −  turbulence model.
Reynolds number a total of 2 074 078 nodes.The angles of the cells are between 30 ∘ and 140 ∘ degrees, with 92% between 70 ∘ and 100 ∘ .The  + values at the walls yield an average value of 50.The continuity equation and the 3D incompressible unsteady Reynolds-averaged Navier-Stokes equations are solved.The eddy-viscosity assumption is used to model the International Journal of Rotating Machinery Experimental

Figure 5 :
Figure 5: Radial and tangential velocities at the impeller outlet at midspan for the steady-state simulation using the  −  turbulence model.

Figure 6 :
Figure 6: Radial velocity in the radial gap (/ 2 = 1.02) at mid-span for the unsteady simulation using the  −  turbulence model.

Figure 7 :
Figure 7: Tangential velocity in the radial gap (/ 2 = 1.02) at mid-span for the unsteady simulation using the  −  turbulence model.

Figure 8 :
Figure 8: Experimental and calculated radial velocities in the radial gap for the unsteady simulation using the  −  turbulence model.

Figure 9 : 3 Figure 10 :
Figure 9: Experimental and calculated tangential velocities in the radial gap for the unsteady simulation using the  −  turbulence model.

Figure 11 :Figure 12 :
Figure 11: Pressure signal for the probe located at the outlet of the impeller (left) and at the outlet of the diffuser (right).The signal in (c) is reconstructed using the 20 main frequencies selected in (b).

Table 1 :
Geometric data and operating conditions.