Analysis of the unsteady flow field in a centrifugal compressor from peak efficiency to near stall with full-annulus simulations

This study concerns a 2.5 pressure ratio centrifugal compressor stage consisting of a splittered unshrouded impeller and a vaned diffuser. The aim of this paper is to investigate the modifications of the flow structure when the operating point moves from peak efficiency to near stall. The investigations are based on the results of unsteady three-dimensional simulations, in a calculation domain comprising all the blade. A detailed analysis is given in the impeller inducer and in the vaned diffuser entry region through time-averaged and unsteady flow field. In the impeller inducer, this study demonstrates that the mass flow reduction from peak efficiency to near stall leads to intensification of the secondary flow effects. The low momentum fluid accumulated near the shroud interacts with the main flow through a shear layer zone. At near stall condition, the interface between the two flow structures becomes unstable leading to vortices development. In the diffuser entry region, by reducing the mass flow, the high incidence angle from the impeller exit induces a separation on the diffuser vane suction side. At near stall operating point, vorticity from the separation is shed into vortex cores which are periodically formed and convected downstream along the suction side.


Introduction
Centrifugal compressors for the aeronautical field are expected to achieve high pressure ratios and high efficiencies at design operating point while minimizing the element size. In this context, typical centrifugal compressor stages are composed of high speed impellers with vaned diffusers to achieve the high pressure recovery in a reduced space. On the other hand, extending the operating range as much as possible is also an important design constraint.
As for axial configurations, the limitation at low mass flow rates comes from the rotating stall and/or surge phenomena. Rotating stall is characterized by the presence of one or several cells rotating around the annulus, either in the impeller or in the diffuser. Surge is a system dependant phenomenon associated to large amplitude oscillations of the pressure through the compressor system [1]. The works of Galindo et al. [2] show that the surge intensity can be modified by changing the length of the duct downstream the compressor. However, operating the system in these unstable conditions induces a dramatic drop of performance associated with mechanical stresses that may cause the failure of the compressor. Therefore, a margin (surge margin) is taken to keep away the operating point from these phenomena, leading to an operating range reduction.
Backswept impellers are widely used to increase the stability of the impeller and finally the stability of the stage [3]. Efforts have also been dedicated to identify flow control 2 International Journal of Rotating Machinery strategies in order to delay the emergence of instability. Different techniques are presented by Skoch [4] in a high speed centrifugal compressor. However, the use of stabilization techniques to extend the operating range requires a good comprehension of flow mechanisms occurring before the stall inception.
Stall phenomenon in compression systems has been studied for almost sixty years. The experimental works of Emmons et al. [5] and Mizuky and Oosawa [6] in a low pressure ratio centrifugal impeller with a vaneless diffuser show that rotating stall occurs in the inducer and leads to surge. In a high pressure ratio centrifugal compressor with vaned diffuser, Wernet et al. [7] and Trébinjac et al. [8] observed separation of the boundary layer on the diffuser vane before the surge inception.
Before the onset of these developed instabilities, two families of precursors have been observed, firstly in axial configuration [9] and then in centrifugal configuration [10]. The first precursor signal is the growth of a small amplitude disturbance with a long length scale referenced as modal stall, while the second concerns the growth of a higher amplitude disturbance but with smaller length scale (several blade passages) termed spike.
Moreover, most of the past works deal either with high pressure ratio (>5) transonic compressors using vaned diffuser or low pressure ratio (<2) compressors using vaneless diffuser. For the first category, there is evidence that the diffuser entry region is often responsible for the surge onset while for the second impeller inducer seems to be the weakest zone in terms of flow instability. The present study concerns a subsonic and moderate pressure ratio (2.5) centrifugal compressor stage composed of a splittered impeller and a vaned diffuser.
In addition, work on the topic comes generally from experimental investigations where the flow description is therefore partial at best. In this study, investigations are performed through a computational fluid analysis. Near the surge line, unsteady phenomena uncorrelated to the blade passing frequency may appear. Therefore, to capture all of the spatial and temporal content, the calculation domain has been extended to all the blade passages leading to 60-millionpoint mesh. Since, the surge onset is generally triggered in the impeller inducer or in the vaned diffuser entry zone, efforts are concentrated at these locations.
The two main objectives of this paper are to (i) numerically analyze the modifications of the flow structure induced by the mass flow reduction from peak efficiency to a near stall condition and (ii) investigate in depth through unsteady flow data sets the near stall operating point in order to propose different scenarios which may lead to flow breakdown. In the first part of the paper, the study case and the numerical procedure are presented. Then, the numerical model is validated thanks to experimental measurements. Afterwards, detailed analysis of the unsteady flow features in the impeller inducer and finally in the diffuser entry zone is given depending on the operating point. The last two parts constitute the scope of the paper.

Test Case Presentation
The centrifugal compressor stage considered for this study has been designed by Liebherr-Aerospace Toulouse SAS and is integrated in an air-conditioning system. The stage is composed of a backswept splittered unshrouded impeller, a vaned diffuser, and a volute. The design specification is based on a stage static-to-total pressure ratio of 2.5 with a design rotation speed of 38000 rpm. The impeller contains 8 main blades and 8 splitter blades with a backsweep angle of 32 ∘ . The impeller exit radius is about 100 mm. The vaned diffuser consists of 21 wedge blades (see Figure 1).

Numerical Procedure
3.1. Flow Solver. Computations are performed using the solver developed by ONERA and CERFACS [11]. It solves the three-dimensional unsteady compressible Reynoldsaveraged Navier-Stokes equations, based on a cell-centered finite volume approach on multiblock structured grids. The turbulent viscosity is computed with the one-equation Spalart-Allmaras model [12]. The convective fluxes are computed with the centered second-order scheme with artificial dissipation of Jameson and diffusive fluxes with a secondorder centered scheme. The time-marching integration is performed with an implicit scheme composed of the backward Euler scheme and a scalar lower-upper symmetric successive overrelaxation method (LUSSOR) proposed by Yoon and Jameson [13]. This time-marching scheme is coupled with a second-order dual time stepping method proposed by Jameson [14]. The number of physical time steps to discretize a complete rotation is set to 1680, corresponding to 210 time steps per impeller main blade passing (impeller has 8 main blades). This discretization value permits a correct description of the impeller-diffuser interactions according to what is usually reported in the literature [15]. For each physical time step, iterations are performed in the inner loop until two orders of magnitude residual reduction are reached. This condition is satisfied in less than 10 subiterations. To reach the periodic state, at least 12 impeller rotations are needed, equivalent to approximately 20000 physical time steps. They are performed using 512 computing cores and require 200000 CPU hours for a single operating point.

Boundary Conditions.
The computational domain contains the centrifugal impeller with the inlet bulb, the vaneless space, and the vaned diffuser and ends with a 90 ∘ turning pipe used as a computational buffer zone to damp possible reflections from the exit boundary condition ( Figure 2). The computational domain extends 1.5 impeller inlet tip radius upstream and downstream of the blade rows. As the meshing of the complete annulus requires significant CPU cost, the volute is not taken into consideration. Considering the inlet boundary conditions, the total pressure, the total temperature, and the flow angles (axial flow) are prescribed. When the operating point moves toward the surge line, the slope of the stage pressure ratio characteristic may reach zero or even positive values. Therefore, a static pressure exit condition is not adapted because two solutions may be International Journal of Rotating Machinery  obtained (different mass flow rates) for a same outlet pressure while a fixed mass flow condition may lead to numerical problem. Therefore, in the present study, the outlet is modeled using a throttle condition, coupled with a simplified radial equilibrium law. The outlet static pressure out is set by the following relation: where 0 is the inlet total pressure,̇( ) is the mass flow rate at iteration through the exit, section and is the throttle parameter. The simulated operating point can move from choked point to surge line by simply increasing the value of the throttle parameter. The rotor-stator interface is treated with a sliding mesh method. In elsA, the communication through the sliding surface is performed using a distribution of fluxes. This approach rigorously ensures conservativity for planar interfaces, which is almost the case here. Details on the implementation and use of the sliding mesh technique with the elsA code can be found in the work of Filola et al. [16] and Gourdain et al. [17]. At the blade, hub and shroud walls no-slip adiabatic conditions are prescribed.

Mesh Parameters.
The structured mesh grid was generated with Autogrid V5 using classical H, O, and C topologies. In order to obtain mesh independent results, the parameters to generate the mesh result from a previous study [18] performed in the same configuration. The size of the first cell is set to 3 m corresponding to a normalized wall distance + well below 3 at the walls. The impeller main blade passage grid and splitter blade passage grid consist of 89 points in the spanwise direction including 29 points in the gap region, 92 points in the pitchwise direction, and 161 points in the streamwise direction. The diffuser blade passage contains 57 points in the spanwise direction, 119 points in the pitchwise direction, and 141 points in the streamwise direction. The impeller blade passage and the diffuser blade passage include, respectively, 2.6 * 10 6 and 1.7 * 10 6 cell grid. The single passage is repeated to obtain the full annulus and the calculation domain reaches a total of approximately 60 million points.

Data Extraction.
Unsteady computations are performed for three operating points OP1 (peak efficiency), OP2, and NS (near stall). As will be illustrated later, unsteady fluctuations for the OP1 and OP2 are only generated by impeller-diffuser interactions. In other terms, for these two operating points, the flow is time periodic in the frame of reference of each row. After reaching the unsteady periodic state, a full rotation of the rotor is performed to extract data (unsteady and timeaveraged). The time-averaging period is equal to one rotor rotation.
For the NS operating point, unsteady effects are not only limited to rotor-stator interactions, and the natural periodicity of the flow is no longer valuable. Therefore, after reaching the stable state, the simulation has been extended during six rotor rotations to validate the stability of the operating point and to extract data. The time-averaging period is therefore equal to 6 rotor rotations.
Because of the domain size, the data extraction of the complete flow field is hardly affordable with a correct temporal resolution. Therefore, data extraction is segregated into local information recorded at each time step, twodimensional planes extracted every 10 time steps, and the complete three-dimensional field saved three times per impeller rotation. Figure 3 depicts the total-to-static pressure ratio defined as = 4 / 01 as a function of the corrected mass flow from numerical simulations results and measurements, for the design speed line. The corrected mass flow is defined aṡ

Numerical Model Validation
The experimental value of 4 is the mean value of three static pressure probes located on the hub surface at diffuser exit radius (plane 4, see Figure 2). The static pressure from numerical results is extracted at the same location. All the computed operating points show good agreement with experimental data. The main objective of the measurements was to validate  the numerical model. In the following sections, investigations will only consider the numerical simulation results.

Impeller Flow Structure Analysis
When the centrifugal impeller is responsible for the surge onset, the origin of the destabilization occurs generally in the inducer and near the shroud. The unsteady numerical results at this specific location are now presented.

Secondary Flow Effects.
The centrifugal impeller flow structure has been studied continually since the 1970s. The work of Eckardt [19,20] shows that the main flow is affected by secondary flows, responsible of the jet-wake structure classically observed in centrifugal impellers. They are produced by reorientation of transverse vorticity (in the boundary layers) into longitudinal vorticity under the effects of curvature and rotation [21]. The meridional curvature part effect drives the low momentum fluid in the boundary layer of the blades (radial migration from hub to shroud), while the rotation effects act mainly on the fluid particles in the hub and shroud boundary layer (migration from PS to SS). The secondary flow effects are generally noticeable from the axialradial turn and intensify up to the impeller exit due to the thickening of the boundary layers. To investigate the impeller flow structure depending on the operating points, analysis of the inlet flow field conditions is first conducted. The operating point displacement from peak efficiency to near stall leads to a gradual decrease of the mean meridional velocity. Since the rotation speed is constant, the incidence angle on the impeller main blades rises. In the study case, the mass flow reduction induces an increase of the incidence angle of approximately 4 ∘ from OP1 to OP2 and 4 ∘ again from OP2 to NS (Figure 4).
Because of the high incidence at low mass flow rate and of the pressure rise occurring in the impeller (adverse pressure gradient), the boundary layer thickens and separates on the impeller main blade suction side. Figure 5 the reduced axial velocity profile in the suction side boundary layer, at 50% of the span near the leading edge for the three operating points. At high mass flow rate (OP1), the incidence angle does not exceed the critical value and no separation is observed. By moving to OP2, the critical incidence angle is reached and a small separation occurs while at NS condition it significantly increases. For OP1, the boundary layer thickness represents 1% of the pitchwise. By moving to OP2, the boundary layer thickness is twice larger while at NS it is four times larger. As a consequence of the boundary layer thickening at low mass flow rate, the secondary flow effects become stronger and are noticeable from the leading edge of the impeller blade. Figure 5(b) plots the reduced radial velocity profile in the suction side boundary layer, at 50% span near the leading edge, for the three operating points. At design condition, the boundary layer is not enough developed (boundary layer is too thin) to induce a radial migration. Therefore, an increase of radial velocity near the suction side of the impeller blade is not observed. At low mass flow rates, due to the thickening of the boundary layer, the secondary flow effects can clearly be observed by the significant rise of the radial velocity near the suction side of the impeller blade. Moving from OP2 to NS leads to intensification of this mechanism and the increase of the radial velocity is larger.
The low momentum fluid near the main blade suction side is then transported along the blade from hub to shroud (positive radial velocity). At the tip of the blade, it is transported and stretched by the leakage flow of the main blade in the middle of the channel ( Figure 6). As a result, by observing the time-averaged meridional velocity at section B (Figure 7), at low mass flow rates, a velocity deficit region (wake) can be noticed in the right channel close to the shroud. This region results from the combination of secondary and leakage flows. By reducing the mass flow along the speed line from OP2 to NS, the wake region significantly enlarges. In the left channel, the mass flow reduction does not affect the flow structure which is principally composed of the core flow.

Unsteady Flow Analysis.
To investigate the unsteady flow pattern for the simulated operating points, a spectral analysis is performed from a local extraction, recorded at each time step leading to an approximate sample frequency of 1 Mhz. In normal operating conditions (stable conditions, constant rotation speed. . .), unsteady phenomena are mainly induced by the blade passing effects. Therefore in the impeller, the flow is time periodic with a period of = 2 /Ω while in the diffuser the period is = 2 /Ω . This specificity is often used to reduce the calculation domain to one single blade passage. Afterwards, the solution can be obtained with an unsteady calculation model using a spatialtemporal periodicity (chorochronic approach). Figure 8 plots the frequency spectra of a static pressure probe linked to the relative frame located at 90% span, at the impeller inlet. It can be seen that the frequency content for OP1 and OP2 is limited to the blade passing frequency located at * = 21 (diffuser has 21 vanes). Since the compressor operates in subsonic conditions, potential effects from the vaned diffuser can propagate upstream and reach the impeller inlet. These operating points could have been simulated using a spatial-temporal periodicity. However, at NS condition, the frequency content is not only limited to the blade passing frequency and a periodic unsteady phenomenon emerges at * = 6. This occurrence shows clearly the need of meshing all the blade passages to the analyzed near stall operating point. The following part focusses on a detail analysis of this frequency emergence occurring at NS.
As discussed in the previous section, the impeller flow structure is composed of the main flow and of the secondary and leakage flows, leading to high and low meridional velocity zones (see Figure 7). The interface between these two flow structures is a region of significant shear. By reducing the mass flow, the meridional velocity deficit due to leakage and secondary flows increases and at NS condition it is such that the velocity gradient is enough to create an interface instability. Vortices are formed at the interface and are transported with the main flow downstream. Figure 9 depicts the instantaneous meridional velocity and instantaneous streamlines at 90% span, in the relative frame. The incoming flow region can be identified as the high velocity zone (blue region) while secondary and leakage flows are marked with a significant low velocity zone (white region). At = 9.1 , a vortex is seen near the splitter blade leading edge while the black point represents the location of the vortex formation. The vortex is then transported by the main flow and grows. At = 10 the vortex is generated and two vortices can be noticed. During one rotor rotation, six vortices are formed and shed, responsible of the fundamental frequency * = 6 seen in the frequency spectra ( Figure 8). All the blade channels create the same phenomenon, but the vortices location varies slightly between the different passages.

Similarities with Axial Compressors.
In axial compressors, März et al. [22] have also observed vortices in the rotor tip region, induced by the interaction between the reverse flow near the trailing edge, the tip clearance flow, and the incoming flow. The vortices move from the suction side to the pressure side and are referenced as rotating instability. Due to the reverse flow near the trailing edge, the vortices do not move downstream and travel circumferentially. The compressor operates in a stable mode even with this rotating instability. In the present case, as there is no reverse flow close to the shroud, the vortices are formed and convected downstream.
According to Duc Vo et al. [23], there are two necessary conditions for the spike disturbance formation in axial compressor. The first one is that the interface between the incoming flow and the tip clearance flow becomes parallel to the leading edge plane, permitting the tip clearance flow to spill into the following blade passage. The second is the initiation of reverse flow at the trailing edge. The first criterion has been investigated in our case by observing the interface position. Figure 10 shows the time-averaged entropy contour map at the tip of the blade. The interface position depends on the flow momentum balance between the incoming flow and the leakage flow. The mass flow reduction induces a decrease of the incoming flow momentum and an increase of the leakage flow momentum due to the blade loading increase. As a consequence, when the mass flow is reduced, the interface between the two flow structures becomes more tangential. For OP1 and OP2, the interface line (red line in Figure 10) goes from the main blade leading edge to the splitter blade leading edge. At NS condition, the interface line is significantly displaced but does not reach the leading edge of the adjacent blade.
Therefore, we hypothesize that if the mass flow is further reduced, the interface between the incoming and leakage flows will become parallel to the leading edge plane, permitting the leakage flow to spill into the following passage and leading to impeller rotating stall.

Diffuser Inlet Flow Structure Analysis
For the present compressor stage, the diffusion process is accentuated through a vaned diffuser. As already mentioned, in the literature concerning configurations with vaned diffusers, there is evidence that surge is often triggered in the vaneless space or in the semivaneless space. This part focuses on the description of the flow at this location.

Mass Flow Reduction Effects.
At the impeller exit, when the operating point is displaced to the surge region along the speed line, the mean tangential velocity increases due to the rise of work input while the mean radial velocity decreases. Therefore, the vanes of the diffuser have to operate with a higher incidence angle. In addition, all along the meridional curvature, secondary and leakage flows continue to interact with the main flow inducing a highly distorted pattern in both the spanwise and the pitchwise directions. The present part investigates the impeller exit flow conditions depending on the three operating points.
Since the work of Dean and Senoo [24], the impeller exit flow structure is described in the pitchwise direction by considering the jet-wake model. The wake zone contains low radial velocity and high absolute tangential velocity flows inducing high absolute flow angle values. Figure 11 profile at 90% span. Considering the three operating points, major differences occur in the right channel (from 50% to 100% of the pitch). The increase of flow angle values at OP2 and NS is linked to the low meridional velocity zone extent at low mass flow rate, noticed in the impeller inducer near the shroud (see Figure 7). Due to the impeller rotation, this pitchwise distortion induces temporal fluctuations of the flow angle at the diffuser inlet, which play a role in the unsteady flow structure described in the following.
Besides this pitchwise distortion, the study of Deniz et al. [25] leads to the conclusion that the diffuser performance is mainly determined by the axisymmetric time-averaged inlet flow angle. Figure 11(b) shows this quantity for the three operating points. Diffuser vane angle (metal angle) is constant from hub to shroud. Due to shroud and hub curvatures, the flow is decelerated and the boundary layer thickness increases on the convex shroud surface although a transfer of flow toward the concave hub side induces a radial velocity increase. Therefore, from hub to 85% span and for the three operating points, the flow angle naturally rises while above 85% span it rudely increases until the shroud due again to secondary flows and leakage flows effects. The compressor stability may be affected when the incidence angle becomes positive ( − Blade > 0), leading to the possibility of boundary layer separation on the diffuser vane suction side. As the mass flow is reduced, the part of the span in that situation significantly extends from 70-100% of the span to 40-100% of the span (see Figure 11(b)).
In addition to the high incidence angle value near the shroud, the pressure gradient in the semivaneless space increases by reducing the mass flow. Therefore, given these conditions, a boundary layer separation occurs on the diffuser vane suction side for OP2 and NS. Figure 12 illustrates the separation zone by representing the time-averaged radial velocity contour and the streamlines. The red zone shows negative radial velocity (reverse flow) while the blue zone shows positive radial velocity. As the flow is subsonic in the vaneless space, the separation bubble yields to a flow deviation, to a rise of flow angle, and finally to a negative radial velocity zone. The mass flow reduction from OP2 to NS induces an enlargement and a displacement of the separation bubble along the suction side toward the leading edge, yielding to a more intense reverse radial flow.

Near Stall Condition Analysis.
The unsteady flow structure in the boundary layer separation shown previously has been analyzed with the 2 vortex criterion [26]. The velocity gradient tensor J is decomposed into its symmetric part S and antisymmetric part Ω and the eigenvalues of J 2 +Ω 2 are determined. If the second eigenvalue 2 is negative in a region, it belongs to a vortex core. This criterion is applied in the vaned diffuser considering the three-dimensional unsteady data. Figure 13 shows the instantaneous negative isosurface of  the second eigenvalue coloured by the normalized span. The separation due the high incidence angle allows the vorticity from the leading edge to be shed to form a vortex-tube. It spans from the diffuser vane suction side (60% span) to the shroud. Due to the large increase of incidence angle from 80% span to the shroud, the upper end (near shroud) of the vortextube tends to move away from the diffuser vane suction side. The recent works of Pullan et al. [27] on axial configuration demonstrate that the process of spike formation is linked to a suction side boundary layer separation resulting from high International Journal of Rotating Machinery incidence angle. A vortex is also formed which starts at the suction surface and terminates at the shroud. It is observed that near the shroud the vortex moves in the circumferential direction along the shroud increasing the flow angle on the adjacent blade. The boundary layer on the adjacent blade separates and the structure can propagate according to the Emmons theory [5]. The onset of spike instability has also been studied numerically by Everitt and Spakovszky [28] in an isolated vaned diffuser. The works show that flow separation at the diffuser vane leading edge associated to radial reverse flow near the shroud allows the vorticity from the leading edge to be convected into the vaneless space. The diffuser inlet blockage rises leading to diffuser instability. This study shows that the high tangential flow at the impeller exit is a key feature for spike onset.
In the study configuration, due to the blade passing effects associated to the high flow angle fluctuations, the vortex core behavior is highly unsteady and time periodic with the blade passing frequency (a blade passage is composed of two channels). Therefore during one channel passing period the vortex forms and expands ( Figure 14). During the second channel passing period (Figure 15), at = 0.0523 the vortex core is detached from the suction side, and the lower end separates from the upper end. As the radial velocity is higher at midspan compared to near the shroud, the lower end of the vortex is convected at higher velocity.
It is found that the vortex-tube does not move in the circumferential direction but downstream along the suction side. Therefore, the scenario described for the spike onset in the literature is not observed for the simulated operating points. However, further mass flow reduction may lead to a more tangential flow (increase of flow angle) near the shroud inducing a vortex displacement toward the circumferential direction and initiating the spike inception process in the vaned diffuser.

Conclusion
Unsteady numerical simulations have been performed in a 2.5 pressure ratio centrifugal compressor stage to achieve a comprehensive description of the flow field from peak efficiency to near stall. The calculation domain extended to all the blade passages has permitted to observe a new unsteady flow pattern with a phenomenon decorrelated to the blade passing frequency when the compressor operates near stall. The conclusions are summarized as follows.
(i) In the impeller inducer, reducing the mass flow induces a rise of the secondary flow effects leading to a significant enlargement of the low momentum fluid region accumulated near the shroud. At near stall condition, the interface between the secondary and leakage flows with the main flow becomes unstable leading to a vortex formation.
(ii) At the blade tip and at near stall conditions, the interface line that demarcates the oncoming flow from the leakage flow is almost aligned with the leading edge plane, suspecting that further mass flow reduction will drive the compressor into stall as observed in axial configurations.
(iii) In the diffuser entry region, the decrease of mass flow rate induces an increase of flow angle leading to a separation on the diffuser vane suction side near the shroud. At near stall condition, the vorticity from the separation is shed and leads to a periodic formation of vortex-tube which travels along the suction side.
(iv) Further mass flow reduction may change the vortex trajectory to a more tangential direction, leading to a propagation of the structure in the circumferential direction and inducing diffuser rotating stall.
(v) Given theses conditions, the prediction of the stall onset region and a priori scenario appear difficult. However, recent investigations presume that the flow breakdown occurs in the impeller inducer due to the alignment of the interface line with the leading edge.