Large-Eddy Simulation of Unsteady Flow in a Mixed-Flow Pump

This article describes the large-eddy simulation (LES) of the internal flows of a high–specific-speed, mixed-flow pump at low flow-rate ratios over which measured head-flow characteristics exhibit weak instability. In order to deal with a moving boundary interface in the flow field, a form of the finite-element method in which overset grids are applied from multiple dynamic frames of reference has been developed. The method is implemented as a parallel program by applying a domain-decomposition programming model. The predicted pump heads reproduce the instability and agree fairly well with their measured equivalents, although the predicted stall takes place at a flow-rate ratio that is a few percentage points lower than the measurements. The phaseaveraged distributions of the meridionaland tangentialvelocity components at the impeller’s inlet and exit cross sections were also compared with those measured by laserDoppler velocimetry. Reasonably good agreements have been obtained between the computed and measured profiles. The developed LES program thus seems to be a promising design tool for a high–specific-speed, mixed-flow pump, particularly for off-design evaluations.

tations of flow fields for engineering applications.Because the RANS equations are stated in terms of time averages, computation with RANS equations has inherent limitations in predicting the unsteady nature of a flow field.Solutions based on the RANS equations usually deteriorate when the flow field of interest involves the large-scale separations that are often encountered in engineering applications.On the other hand, the large-eddy simulation (LES), in which turbulent eddies of a scale larger than the computational grid are directly computed, has the potential to predict unsteady flows and flow fields that include regions of large-scale separation much more accurately than RANS-based computation does in general.The ultimate goal of the present study is to develop a practical engineering tool that is based on LES, with a particular emphasis on internal-flow simulations of turbomachinery and simulations of aeroacoustics (Kato et al., 1991(Kato et al., , 1993(Kato et al., , 1999(Kato et al., , 2000)).
For turbomachinery simulations, it is necessary to deal with moving boundary interfaces between flow fields, such as those that appear between a rotating impeller and a stationary casing.In the present study, the moving boundary interface is taken into account by using overset grids from dual frames of reference.In this method, the overall grid is composed of several (usually between two and five) grid sets, and appropriate transactions take place at the interface regions.In consideration of its applications to a complex geometry, the finite-element method is used to discretize the governing equations of the flow field.It is implemented for parallel processing, and it is therefore possible to complete a large-scale flow-field computation within a practical period of computation on a current-model distributed-memory parallel computer.
For simulations of unsteady flows in turbomachinery, this method was first applied to the prediction of internal flows in a mixed-flow pump (Kato et al., 1999).The predicted fluctuating forces on the impeller as well as the static pressure in the casing agreed quantitatively with their measured equivalents.This article describes applications of this method to the internal flow of a mixed-flow pump stage with a high specific speed.A particular emphasis is placed on the capability of this method to predict unstable head-flow characteristics.Phase-averaged distributions of the impeller's inlet and exit velocities are compared with the values measured by a two-dimensional laser-Doppler velocimetry (LDV).

GOVERNING EQUATIONS
The governing equations on which the present study is based are the spatially filtered continuity equation (Equation ( 1)) and the Navier-Stokes equations (Equation ( 2)) for the flow of an incompressible fluid, as represented in Cartesian coordinates: where ūi (i = 1, 2, 3) is the grid-scale velocity component in the x i direction, p is the grid-scale static pressure, ρ is the density, and ν is the kinematic viscosity.f i is the inertial force associated with the motion of the frames of reference.In particular, for a stationary frame of reference, whereas for a rotational frame of reference, the centrifugal forces and Coriolis forces must be added: where is the angular velocity of the frame of reference, which is assumed to be rotating about the positive x 3 axis.The effects of those eddies that are not resolved by the grid (subgrid-scale eddies) are modeled after Smagorinsky (1963), incorporated with the Van Driest wall-damping function that represents the near-wall effects: The model coefficient C S is fixed to 0.15, which is a standard value for flows with large separation, and the grid-filter size is computed as the cube root of the volume of each finite element.

NUMERICAL METHOD Overset Grids from Multiple Dynamic Frames of Reference
In the present study, a moving boundary interface in the flow field was treated with overset grids from multiple dynamic frames of reference (Kato et al., 1999).The application of this method to the interaction between a rotating impeller and a stationary casing in a pump is schematically depicted in Figure 1.A computational mesh that rotates along with the impeller is used to compute the flow within the impeller.On the other hand, a dedicated stationary computational mesh for each part computes the flow in stationary parts of the pump.Each mesh includes appropriate margins of overlap with its neighboring meshes upstream and downstream.At every time step, the velocity components and static pressure within each such margin are the values interpolated in the computational mesh of the corresponding neighbor.Elementwise trilinear functions are used to interpolate both the velocity components and the static pressure (Ikegawa et al., 1994;Kaiho et al., 1997).When velocity components are overset, an appropriate coordinate transformation must be applied to take into account the differences between the frames of reference.

Finite-Element Formulation
A streamline-upwind finite-element formulation (Kato et al., 1991) is used to discretize the governing equations of the flow field.This formulation is based on the explicit Euler's method, but it shifts the spatial residuals of the governing equations in the upstream direction of the local flow.The magnitude of this shift is one half of the time increment multiplied by the magnitude of the local flow velocity.This shift exactly cancels out the negative numerical dissipation that is otherwise the result of applying Euler's method, and it guarantees stability and the accuracy of the solutions.The proposed formulation essentially possesses second-order accuracy in terms of both time and space

FIGURE 2
Sustained parallel computing performance.and has been successfully applied to the LES of external as well as internal flows (Kato et al., 1999(Kato et al., , 2001)).

Parallel Implementation
The global computational domain is partitioned into a prescribed number of subdomains, and each of the subdomains is assigned to a dedicated processing node.The communicating pairs and local coordinates where given interpolations are taking place change at the moving-boundary interfaces at every time step.The usual (unsophisticated) parallel implementation therefore includes broadcast communications at every time step as the processing node searches for its new communication pairs.This communication overhead seriously deteriorates the overall parallel computing performance.In the present study, the communication pairs were searched for in advance by a serial computation and were fed to the parallel flow solver as input data at every time step.This procedure not only avoids the otherwise inevitable communications overhead but also brings in greater flexibility with the flow solver.Figure 2 is a plot of the overall performance sustained on Hitachi's SR8000 super-computer against the number of the processing nodes.A parallel computing efficiency of over 85% is achieved on this platform.

SIMULATION RESULTS
The method was applied to the prediction of the internal flows of a mixed-flow pump stage that are composed of an inlet whirl stop, a four-bladed, open-shrouded mixed-flow impeller, and a discharge casing with eight diffuser vanes.This pump has a designed specific speed ω S of approximately 3.0.The measured head-flow characteristics of this pump exhibit weak instability over a flow-rate ratio of 50-60%.The authors' interests were in whether the LES could predict this instability.The phaseaveraged distributions of the meridional and tangential velocities at the impeller's inlet and exit cross sections were measured by a two-dimensional LDV, and the results were also compared with the computations.

FIGURE 3
Computational mesh for a mixed-flow pump stage, composed of meshes for five parts.Note: The mesh for the inlet driving section is not shown here.
The computational mesh used in this study is shown in Figure 3. Upstream of the regulation-plate mesh, a driving section (not shown in this figure) is overset and therefore, at every time step, the inlet velocities of the regulation-plate mesh were the values interpolated from instantaneous velocities of the driving section (see Figure 4 for fluctuations in centerline velocities in this section).The flow rate of the pump was set by adjusting the magnitude of uniform acceleration that was applied to the driving section in the axial direction.The Reynolds number based on the impeller's exit diameter D 2 and the circumferential velocity u 2 at this point was 5.7 × 10 6 .In order to investigate the

FIGURE 4
Fluctuations in centerline velocities in the driving section.

FIGURE 5
Comparison of measured and predicted head-flow characteristics.
effects of grid resolution on the accuracy of the solution, two sets of meshes with different grid resolution were used.The numbers of elements were approximately 1.7 million for the coarse mesh and 5 million for the fine mesh.The actual impeller possesses a tip clearance of approximately 1.0 × 10 3 D 2 , but this clearance was not taken into an account in the present study because of the limitation of available grid number.The time increment was set so that 8000 time steps corresponded to a single revolution of the impeller and the total pump head and mean-velocity distributions were evaluated by averaging the flow field during 10 revolutions of the impeller.
The computed total pump-heads, together with their measured equivalents, are plotted in Figure 5.The total pump-heads predicted by the fine-mesh LES agree fairly well with the measured values down to a flow-rate ratio of 60%, where the measured head-flow characteristics indicate an onset of the total stall.The fine-mesh LES predicted this onset at a flow-rate ratio of 54%, which is 6% lower than the measured value.The maximum difference in the predicted and measured total heads was about 5%.The discrepancy between the computed and measured characteristics may be attributed to (1) the grid resolution around the blade's surfaces being insufficient to resolve the turbulent boundary layer developing on them; (2) the effects of the leakage flow through the tip clearance that were not considered in the present simulation; and (3) the subgrid scale model adopted in the present study that was not capable of predicting transition to turbulence.Further investigations to identify the cause of this discrepancy are currently under way.
Here, we present flow fields at flow-rate ratios of 60% and 43% because they respectively represent the pre-stall and poststall flow fields.Figure 6 shows the typical instantaneous distributions of the axial-velocity component at the impeller's inlet cross section as well as at the center plane of the pump.In this figure, regions that are colored magenta indicate recirculation of the flow.At the flow-rate ratio of 60%, small recirculations are already taking place near the leading edges of the impeller's blades, around their tip.As the flow-rate ratio is reduced to 43%,

FIGURE 6
Instantaneous distributions of axial velocity at the impeller's inlet cross section and center plane.(a the above-mentioned recirculations are considerably enhanced and their regions are extended in both pitchwise and spanwise directions.Note, however, that the instantaneous recirculations are still local, as concentrated around the leading edge of the impeller's blades.Local recirculations are also taking place at the inlet of the diffuser, around its hub side, at both flow-rate ratios, but they are turned to the centrifugal direction at the lower flow-rate ratio.
Figure 7 shows instantaneous distributions of the tangential velocity component of the absolute flow at the same instance and cross section as is shown in the previous figure.In this figure, regions colored blue indicate that the local flow is rotating in the direction identical to that of the impeller's revolution, while regions colored orange are rotating in the opposite direction.Small regions near the leading edge of the blades around their tip are rotating in the identical direction to the impeller's revolution even at the pre-stall flow-rate ratio.Further examination of the time evolution of the instantaneous flow field revealed that these regions are rotating at about 60% of the impeller's tip speed, thus propagating, when they are viewed from the rotational frame of reference, in the backward direction at a speed roughly corresponding to 40% of the tip speed.This rotation is confined to a small portion near the tip, but its origin seems to be similar to that of the rotating stall seen in axial compressors.
The vortex structure near the tip is now being investigated in detail.
When the flow-rate ratio is reduced to 43%, a substantial portion of the flow near the shroud starts to rotate in the same direction as the impeller's revolution.This large region of prerotation apparently leads to reduction in the total pump head, which results in instability in the head-flow characteristics.
As previously mentioned, the phase-averaged distributions of the meridional-and tangential-velocity components were measured by a two-dimensional LDV at the impeller's inlet and exit cross section.Figure 8 compares the measured (left) and predicted (right) meridional velocity at the impeller's inlet cross section for the pre-stall (60%) and post-stall (43%) cases.The

FIGURE 8
Comparisons of phase-averaged meridional velocity at the impeller's inlet cross section (left, LDV; right, LES).
overall agreement between the computations and measurements seems fairly good for both flow-rate ratios.
Figure 9 shows comparisons similar to those in the previous figure, where the tangential-velocity component of the absolute flow is compared.Again, the predicted distributions are in good agreement with their measured equivalents for both flow-rate ratios.The small region near the shroud that has the forward rotation seen in the instantaneous flow field at the pre-stall condition can still be identified in the phase-averaged flow field.The measured distributions are somewhat more uniform in the pitchwise direction than in the computations.This discrepancy between the computations and the measurements may be attributed partially to the difference in the method applied to taking a phase average.The computed result is of time-average in the rotational frame of reference, thus averages flow at a position in the impeller, over all the impeller's angles, that is relative to the casing.On the other hand, the measurement averages flow at a position in the impeller that is based on its particular angle relative to the casing.Figure 10 compares predicted and measured pitchwise average distributions of the meridional-and tangential-velocity components at the impeller's inlet cross section, where C m and C u , respectively, denote the meridional-and tangential-velocity components.The velocity profiles predicted by the fine-mesh LES are in good agreement with the measured profiles except near the shroud.For the pre-stall case (60%), the computed profile of the tangential-velocity component shows weaker prerotation near the shroud than the measurements show; this is the result of the smaller amount of recirculation predicted by the computation.This discrepancy between the computations and the measurements is directly related to the difference in the flow-rate ratio at which the total stall takes place.The large portion of flow, however, that is recirculating and that results in the intense prerotation, is accurately predicted by the present LES.Because the essential feature of LES is to directly compute turbulent eddies of a scale larger than the computational grid,

FIGURE 12
Comparisons of phase-averaged distributions of meridional velocity at the impeller's inlet cross section for Q/Q d = 43%.the accuracy of LES becomes, in general, better when the active eddies in the flow become greater in size.In that sense, we can expect a more accurate solution of LES in the post-stall condition (than in the pre-stall condition), where large-scale separation takes place at the impeller's inlet.Finally, Figure 12 shows the computed and measured phase-averaged distributions of the meridional-velocity component at the impeller's inlet cross section in the post-stall case.The predicted profiles are in fairly good agreement with their measured equivalents for all the spanwise positions in which the comparison was made.

CONCLUSIONS
Large-eddy simulation was applied to the prediction of internal flows in a high-specific-speed, mixed-flow pump stage that possesses weak instability in its head-flow characteristics at low flow-rate ratios.The standard Smagorinsky model, incorporated with the Van Driest damping function, was adopted as the subgrid-scale model.The numerical method was based on a streamline-upwind, finite-element formulation with accuracy of the second order in both time and space, and it incorporates the application of overset grids from multiple and dynamic frames of reference.The method is implemented as a parallel program by a domain-decomposition programming model and, therefore, it can be applied to a large-scale computation on a distributedmemory parallel computer.
The predicted head-flow characteristics were in fairly good agreement with their measured equivalents, although the LES predicted the stall point at a somewhat lower flow-rate ratio than did the measurements.The phase-averaged distributions of the meridional-and tangential-velocity components were also compared with those measured by an LDV.Reasonably good agreements were obtained between the computation and the measurement.
The proposed method thus seems to be a promising candidate for use as a hydraulic design tool in the next decade, particularly for off-design evaluations.

FIGURE 1
FIGURE 1 Schematic view of an example of overset grids from dual frames of reference.