Abstract
The harmonic balance method was used for the flow simulation in a centrifugal pump. Independence studies have been done to choose proper number of harmonic modes and inlet eddy viscosity ratio value. The results from harmonic balance method show good agreements with PIV experiments and unsteady calculation results (which is based on the dual time stepping method) for the predicted head and the phase-averaged velocity. A detailed analysis of the flow fields at different flow rates shows that the flow rate has an evident influence on the flow fields. At 0.6, some vortices begin to appear in the impeller, and at 0.4 some vortices have blocked the flow passage. The flow fields at different positions at 0.6 and 0.4 show how the complicated flow phenomena are forming, developing, and even disappearing. The harmonic balance method can be used for the flow simulation in pumps, showing the same accuracy as unsteady methods, but is considerably faster.
1. Introduction
Flow in centrifugal pumps produces a complex three-dimensional phenomenon involving turbulence, secondary flows, separations, and so forth (Brennen [1]). For a low specific speed centrifugal pump, the geometry is complex and asymmetric due to the volute shape; the relative movement between impeller and volute generates an unsteady interaction which not only affects the overall pump performance (flow structure, losses) but also is responsible for pressure fluctuations.
A number of researchers have contributed to the understanding of complex unsteady phenomena by using CFD [2–5]. Although predictions on the unsteady flow field are always valuable, numerical simulation of centrifugal pumps is not easy due to the usual CFD difficulties: turbulence modeling, flow separation, boundary layer, and so on [6]. Besides that, there are also specific problems, as the following:(i)Complex geometry: structured mesh is not easy to generate and difficult to converge; a great number of cells are needed to get more details on the curved blades.(ii)The interaction between impeller and volute requires an unsteady solution process to calculate the time dependent flow field. In addition, the blade position with respect to the volute tongue must be taken into account. This can be partially accomplished in a quasi-steady way: calculating the steady solution with different grid positions. Nonetheless, it is always much better if the code is able to perform an unsteady flow calculation at the same time that slides the impeller grid at each time step.
As computer power increases, increasingly sophisticated computer fluid dynamic models of unsteady flows about turbo machines have been developed. Some researchers have performed time-accurate Large Eddy Simulations (LES) of unsteady flows about pumps [7]. This method provides more accurate solutions than uRANS. However, such codes are expensive to use, requiring significant computer resources. This will make them difficult to use in the design phase.
To analyze the nonlinear hydrodynamics of pump impeller, in the turbo machinery community, the need for efficient unsteady flow solvers led investigators to develop efficient mixed time and frequency domain techniques. Hall et al. [8–10] proposed the harmonic balance method. This method can be used if the flow field is periodic in time, a condition satisfied to a good approximation for many unsteady flows of interest. The unsteady flow field is transformed into a coupled set of steady-state equations. Because time does not appear explicitly, the harmonic balance equations can be solved very efficiently using the same numerical algorithms developed for steady-state flow problems using convergence acceleration techniques.
A number of investigators have used this technique to solve interesting fluid dynamic problems. Such problems include the vortex shedding of a circular cylinder [11, 12], multistage turbo machinery analysis [13, 14], limit-cycle oscillations of wings [15], or normal synthetic jet in quiescent background flow [16]. Recently, investigators [17, 18] have started using similar techniques to analyze unsteady aerodynamics of helicopter rotors and have demonstrated that accurate solutions can be obtained efficiently. In this work, we apply the harmonic balance method (HBM) to a low specific centrifugal pump to demonstrate that these unsteady flows can be accurately calculated with substantially less computational effort than conventional time-accurate solution methods. To the authors best knowledge it appears to be the first time that the HB method has been used to a radial pump.
1.1. The Harmonic Balance (HB) Method and Computational Scheme
Consider the three-dimensional unsteady Navier-Stokes equations. In Cartesian coordinates, these equations are given bywhere is the vector of conservation variables. The impeller rotates with a frequency . Therefore, the flow is periodic in time, with period . Because the flow is temporally periodic, the flow variables may be represented by a Fourier series in time with spatially varying coefficients. For example, the conservation variables may be expressed as an infinite Fourier series given bywhere , , and are the Fourier coefficients of the conserved flow variables.
The Fourier series is truncated to a small number of harmonics :
We denote this vector of flow variables by . Therefore, (3) can be written in matrix form aswhere is the inverse discrete Fourier transform operator.
The Fourier series is then substituted back into the Navier-Stokes equations and rearranged into the individual harmonics. Consider
This gives a system of equations for the Fourier series coefficients. Since the only way for the equations to be true is that every equation set must be zero:
Using this transformation (discrete inverse Fourier transformation) and substituting this into the equations in the frequency domain give
This gives a system of equations for the conservative variables:
is a spectral approximation of the time derivative and is very accurate. In order to solve this system of conventional steady-state equations one integrates in pseudotime :where is a fictitious or pseudotime, used only to march (9) to steady state, driving the pseudotime term to zero. Note that pseudotime harmonic balance equations are similar in form to the original time domain form of the Navier-Stokes equations. Thus, existing well-developed steady CFD techniques may be used to efficiently solve the nonlinear harmonic balance equations, with a comparable number of iterations required. Also, because only steady-state solutions are desired, local time stepping and multiple grid acceleration techniques are used to speed up convergence.
The method has some similarities to the dual time step method, used by Davis et al. [19], Sayma et al. [20], and others to compute unsteady flows in the time domain. This approach has a number of important differences. First, in the dual time step method, one marches from one time level to the next time accurately, using pseudotime steps to drive the residual of the time-accurate equations to zero. The process is repeated over many time steps for several periods until a periodic solution has been reached. Second, because we solve for the solution over one complete period, a spectral operator may be used to compute the physical time derivative . The spectral time derivative is much more accurate than finite difference operators, which are used in the dual time step approach. Therefore, many fewer physical time levels are required using the present method.
The harmonic balance method is implemented into the in-house code SPARC [21]. The equations used are those for compressible fluids. Since we want to use this approach later on for the simulation of compressible cavitation, we have to prescribe the liquid by a suitable equation of state. In this application, we have used the stiffened gas EOS (Tamman), which relates the pressure, density, and temperature by [21]where and are empirically determined constants for the used liquid.
In order to solve for these numerically very stiff equations we are using the preconditioning method proposed by Shouqi et al. [22]. Since the inclusion of the spectral approximation of the time derivative adds a further eigenvalue into this coupled NS system, we have taken that into account by adding a modification according to van der Weide et al. [23]. The numerical scheme used is a 2nd order accurate central scheme with artificial dissipation (switch scheme) according to Jameson et al. [24]. A three-stage Runge-Kutta method with implicit residual averaging, local time stepping, and multigrid acceleration is used to overcome the stiffness of the NS equations.
1.2. Model Description and Computational Setup
The shrouded impeller having backward curved blades with a vaneless single tongue volute obtained from the optimization method described in [25] has been scaled down from a full scale to a model scale. In this scaling, the flow coefficient and head coefficient are retained, whereas the Reynolds number in the impeller is reduced by a factor of 2 by reducing the rotating speedy from 2900 rpm to 1450 rpm. All the geometry parameters are retained, as shown in Table 1. The transparent impeller and the volute with rectangular outside section have been manufactured in Perspex (polymethyl methacrylate (PMMA)) and have been installed in an open test rig, as shown in Figure 2.
1.2.1. Geometry Modeling
Due to the asymmetric volute, the flow in one blade flow passage is not periodic, so the computational domain includes a rotational zone (the whole impeller and the flow of inlet tube, while the wall of the inlet tube should be set as fixed wall) and stationary zones (volute and outlet tube). For numerical stability reasons, and to minimize boundary condition effects, this computational domain is extended upstream and downstream. The inlet tube with the length of five times the diameter is also transparent, and the outlet tube is extended to the length about three times the diameter of the outlet tube, which are corresponding with the position of pressure sensors, respectively. The CAD of this configuration, as well as the boundary conditions, is shown in Figure 3.
1.2.2. Mesh Generation
Structured H–O–H grid topology cells are generated to define the impeller, the inlet tube, and the volute (together with the outlet tube). The geometry has been simplified in the front wearing ring; that means that the leakage flow through it is not considered. A 3-level full multigrid technique was used, the mesh refinement zones are defined near the blades inlet and the volute tongue, and the grids numbers are shown in Table 3. A view of the generated grid at the 2nd finest level can be seen in Figure 4.
1.2.3. Boundary Conditions
Mass flow rate condition at the inlet and static pressure outlet boundary conditions has been used, according to the values measured in the experiment (Figure 3). A sliding mesh technique is applied to the interfaces (impeller-volute interfaces in Figure 3) in order to allow the unsteady interactions between the impeller and the volute. On the surface of the impeller and the volute, no-slip and isothermal wall boundary conditions are applied. In the case of the impeller the wall velocity has been set according to the rotational frequency and a roughness value of 20 μm has been used in order to take the extra losses into account.
1.2.4. Numerical Setting
The calculations were done using our in-house developed code called SPARC. The NS equations were solved, respectively, with the HB method and the dual time stepping method. We have been using 2 harmonic modes giving time slices, as shown in Figure 1 (points 1–5), for the simulation.
The turbulence was modeled with the Spalart-Allmaras one-equation model. The harmonic balance calculations have been done with a three-stage Runge-Kutta scheme (CFL number of about 2.2) using full multigrid, implicit residual smoothing, and local time stepping to accelerate the convergence rate considerably. For the unsteady calculation, a second-order accurate dual time stepping approach was used. All convergence acceleration techniques used for the HB method were also applied as well. The physical time step used here was s, this means that we have done 4000 time steps for resolving one revolution. The number of time step is larger than the usually taken 300–400 time steps, to accurately resolve the physical effects in time.
2. Independence Analysis
2.1. The Selection of the Harmonic Modes
2.1.1. Predicted Head for Different HB Modes
To select an economic harmonics value is considered next. 1, 2, and 3 harmonic modes were investigated based on the 2nd level mesh. The calculated mean head with different harmonic modes and their CPU time is shown in Table 2. It is compared with the value from unsteady calculation (based on the dual time stepping method). From the results, we can extract that for increasing number of harmonic modes, the calculated heads are converging toward the value from the unsteady simulation. The largest error is, as expected, when using one harmonic mode and shows a relative error of 3% only. Not surprisingly, the dual time step solver is the most computationally expensive. It needs about 20 times more cpu-time compared with HB using three modes. Since the predicted head varies a little when we compare the 2 modes with the 3 modes results for this flow rate, as well as for another flow rate (not shown here), we concluded that 2 modes are a good compromise for the further investigations.
2.1.2. Velocity Comparison with PIV
To verify the results from HB method, we compared the absolute velocity inside the volute with the results from unsteady calculation and PIV. The absolute velocity is compared at the same radius ( mm) along an arc between 20° and 120° at the midheight of the impeller outlet as shown in Figure 5(a). The absolute velocity results from PIV and HB are ensemble averaged, while the results from the dual time stepping are time averaged, as shown in Figure 5(b).
(a) The arc line
(b) The mean absolute velocity results
The dimensionless velocity is approximately 0.5 in a large part of the channel. Close to the volute tongue (angle about 20°), there is a zone with low absolute velocity. The results from the PIV test are uniform compared with those from CFD. With just one mode we can get similar results with PIV, while with the increasing of the number of modes, the results are converging toward that of the unsteady calculation.
The difference between the unsteady calculation and the harmonic balance method is largest in the range of the angle between 80° and 120°. There the unsteady calculation shows a larger dimensionless velocity than the harmonic balance one. Since the mesh, the numerical scheme, and the turbulence model used are the same, only the number of modes used for the comparison makes the difference. In the case of a pump with a spiral volute, we must calculate the whole impeller and volute. That means that the time slices used by the harmonic balance have to be equally distributed around 360°. When averaging these results for the harmonic balance method we have only a small number of solutions available compared with a large number in the case of dual time stepping method. Therefore, we believe that increasing the number of harmonics will give a better prediction for the velocity distribution.
2.2. Mesh Independence
2.2.1. Value and Refined Mesh
From the calculations, we got the values distribution of the first mesh. The values at the blades’ pressure side were between 20 and 50 and the values at the suction side were around 20, while the values for all the walls were found to be bigger than 100. Because we did not use wall functions, we need to get proper resolution of the boundary layer. The meshes at the boundary have to be refined by about 10 times.
The number of grid points of the refined mesh is shown in Table 3. The distribution around the blades for the new mesh at the 1st level is about 10 and those at the 2nd level around the blades are less than 5 and for the finest mesh, is in the order of 1.
The number of the cells inside the boundary layer is sufficient according to our experience, it is confirmed by the predicted results for the pump performance, and so it allows the analysis of the main flow phenomena involved.
2.2.2. Mesh Independence Analysis
One mode was chosen to verify the mesh independence, as shown in Figure 6. Please note that the inverse logarithmic scale of the mesh number is shown in this figure. From the figure, we can see that the head value increases with increasing mesh numbers. The coarsest mesh shows rather large disagreement compared to the experiment, while the predicted head from the second level and the third level changes slightly and is close to the measured one. Although the used mesh does not show perfect mesh independence, we have used the second mesh as the basis for the comparison between fully unsteady calculation and harmonic balance method. For a detailed comparison with the experiments, we will use the finest mesh in a forthcoming paper.
2.3. Specification of the Inlet Turbulence Quantities
The ratio is the ratio between the turbulent viscosity and the molecular dynamic viscosity . Unfortunately, no turbulence quantities have been measured at the experiment. In order to prescribe a realistic inlet boundary condition for the turbulence variables in the Spalart/Allmaras model, we need to estimate the eddy viscosity ratio. We have investigated the influence of the results by varying the value of the turbulence viscosity ratio from 5 to 100 for the case with 2 HB modes, and, for the second level mesh, the head changes just about 0.083%. Therefore, in the following calculation, the turbulence viscosity ratio at the inlet was set as 100.
3. Results Analysis
3.1. Predicted Performance
Figure 7 shows the predicted head as a function of the dimensionless flow rate compared to the test results. At flow rates close to the design flow rate and above it (0.8–1.4), the predicted head agrees very well with the tests. At very low flow rates (<0.4) it is well known that the predictive accuracy depend mainly on the turbulence model used. The reason for it is that the flow field becomes more and more rotational in combinations with strong curvature effects and linear eddy viscosity models having difficulties in predicting this flow fields. At flow rates 0.6 and 0.8 the calculated head is somewhat smaller than those measured in the experiment. The reason for it is not clear since little separation effects prevail there. But it is certainly not a problem of the harmonic balance method since an unsteady calculation with dual time stepping shows a similar behavior. We would expect an improvement of the accuracy when using the Spalart-Allmaras turbulence model with a rotation/curvature correction method as proposed by Zhang and Yang [26].
3.2. Flow Analysis at Different Flow Rates
3.2.1. The Velocity at the Design Point and Comparison with PIV Results
In Figures 8(a) and 8(b) the dimensionless absolute velocities at and are compared with the PIV results, respectively. Inside the impeller and volute, the velocities from CFD and PIV have similar value and distribution laws. For the absolute velocity at , the averaged velocity calculated from HB method simulation is very close to the result from unsteady calculation and PIV, except the regions near the volute tongue.
(a)
(b)
3.2.2. Absolute Velocity Contour at Different Flow Rates
The contours of absolute velocities at different flow rates for the first position are shown in Figure 9. For lower flow rates such as 0.4 and 0.0, the flows mainly rotate inside the impeller with very high velocity. In addition, at the zones of the impeller outlet, there are the maximum absolute velocities, with large unsteady flow characteristics. With increasing flow rates, the distributions become more uniform, and an increasing flow rate is convected to the outlet tube, which also contribute and add some dynamic energy to the head. Meanwhile, from the analysis prescribed above, we know that the absolute velocities inside the impeller and volute coincide with experiments very well, except the zones around the impeller-volute interface, which are somewhat higher compared to test results.
3.2.3. Contour of Relative Velocity at Different Flow Rates
Figure 10 shows the contour of relative velocity at different flow rates. At the design point, the flow is not separated and almost uniform. At the flow rate of 0.8, some small vortexes begin to appear at the blades pressure side around the trailing edge just in two or three flow passages. When the flow rate is further decreased to 0.6 the vortex strength increase, especially in the flow passage closest to the volute tongue. At the flow rate of about 0.4, the impeller is almost blocked by six vortexes, and additionally some small vortexes appear inside the impeller. On the other side, at the larger flow rate 1.4, the flow is uniform in most flow passages, except the one opposite to the volute tongue; there is a small vortex visible near the blade inlet.
3.2.4. Contour of Eddy Viscosity Ratio at Different Flow Rates
Eddy viscosity ratio is shown in Figure 11, it is a quantity indicating the strength of the turbulence. For the smaller flow rates (from 0.0 to 0.6), corresponding to the location of vortices in Figure 10, there are some zones with high eddy viscosity inside the impeller, and there is a zone with maximum eddy viscosity in the volute opposite to volute tongue. There are also some cores with high eddy viscosity inside the volute, which are located around each blades trailing edge. Increasing the flow rates, the values of the eddy viscosity inside the impeller become lower and the cores of the eddy viscosity inside the volute shrink and move out to the outlet pipe gradually.
3.3. Inner Flow Analysis at Different Positions
3.3.1. Absolute Velocity Distribution at 0.6
Harmonic balance method provides flow information at positions for each simulation. By analyzing the relationship between the absolute velocities in Figure 9 and the relative “streamlines” in Figure 10, we can deduce that the relative vortexes always appear at the area with high absolute velocity value in the range of 11-12 m/s. The absolute velocity distribution at the other four positions is shown in Figure 12. The main differences appear around the area of the impeller-volute interface. In general, the area with higher velocity exists at the corner of blades trailing edge suction side and due to the asymmetric volute at the corner between blade pressure side and tongue we observe an area with higher velocity.
3.3.2. Relative Velocity
Figure 13 shows the relative velocity and some visualization of it at the other four positions at 0.6. and 0.4, respectively. Corresponding to the area with higher velocities in Figure 12, especially at the blades pressure side, several vortices exist or are forming and disappearing at 0.6, and the flow at the interface begin to become more unstable. Obviously, the starting point for the creation of the separated flow is at the tongue. Here the flow accelerates at the tongue and behind the tongue the flow decelerates and separates. For the 0.4, there are six vortices blocking the flow passage and some small vortices appearing inside the impeller, which contributes to the asymmetrical flow characteristics.
(a) 0.6
(b) 0.4
4. Conclusions
In this work, harmonic balance method has been applied to a low specific speed radial volute pump to predict the unsteady flow fields. The following conclusions can be drawn.(1)Some works were done investigating mode independence of the harmonic balance method. The second finest mesh and suitable turbulence boundary conditions have been chosen for the investigation of the performance curve. It resamples a good compromise between computational effort and accuracy.(2)The CFD results show good agreements with PIV experiments for the predicted head.(3)A detailed analysis of the flow fields at different flow rates shows that the flow rate has a very evident influence on the flow fields. At 0.6, some vortices begin to appear in the impeller, and from 0.4 to 0.0 some vortices block the flow passage.(4)The flow fields at different positions at 0.6 and 0.4 show how the complicated flow phenomena are forming, developing, and even disappearing.(5)The harmonic balance method predicts very similar flow fields like an unsteady simulation but being more than ten times faster.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research has been done while the second author visited the ISTM of KIT, Germany. She has been supported by the National Natural Science Foundation of China (no. 51009072) and the Foundation for senior person with ability in Jiangsu University of China (no. 08JDG040), and a Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions.