Next Article in Journal
Parameterised Model of 2D Combustor Exit Flow Conditions for High-Pressure Turbine Simulations
Previous Article in Journal
Investigation on the Flow in a Rotor-Stator Cavity with Centripetal Through-Flow
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Analysis of Vaneless Diffuser Stall Instability in a Centrifugal Compressor †

1
Royal Institute of Technology (KTH), Department of Mechanics, Competence Center for Gas Exchange (CCGEx), Osquars backe 18, SE-10044 Stockholm, Sweden
2
Baker Hughes, a GE company, 50127 Florence, Italy
*
Author to whom correspondence should be addressed.
This paper is an extended version of our paper published in Proceedings of the European Turbomachinery Conference ETC12 2017, Paper No. 175.
Int. J. Turbomach. Propuls. Power 2017, 2(4), 19; https://doi.org/10.3390/ijtpp2040019
Submission received: 3 October 2017 / Revised: 21 November 2017 / Accepted: 24 November 2017 / Published: 28 November 2017

Abstract

:
Numerical simulations based on the large eddy simulation approach were conducted with the aim to explore vaneless diffuser rotating stall instability in a centrifugal compressor. The effect of the impeller blade passage was included as an inlet boundary condition with sufficiently low flow angle relative to the tangent to provoke the instability and cause circulation in the diffuser core flow. Flow quantities, velocity and pressure, were extracted to accumulate statistics for calculating mean velocity and mean Reynolds stresses in the wall-to-wall direction. The paper focuses on the assessment of the complex response of the system to the velocity perturbations imposed, the resulting pressure gradient and flow curvature effects.

1. Introduction

Ideally, for an efficient design the centrifugal compressor provides a steady fluid flow at an elevated pressure with minimal losses. However, the performance of the centrifugal compressor is harmed under the off-design operating conditions occurring at low flow rates. In such circumstances flow instabilities are developed through the compressor, which may cause large fluctuations in the flow variables and even introduce vibrations with high structural stress levels into the machine. This is an unwanted situation and the aim is therefore to understand the physical mechanisms responsible for the developed flow instabilities. Through linear stability analysis, see e.g., [1], the diffuser rotating stall is said to occur when the local return flow angle relative the meridional exceeds 90 . However, a diffuser instability generally develops earlier than 90 , i.e., prior to purely tangential flow. There exist some correlations that indicate the instability angle to be around 75 –85 with respect to the radial direction depending on the operating conditions, see e.g., [2]. The theory [1] also states that it occurs with steady, uniform inflow and outflow conditions and thus independent of the interaction with the rotating impeller. A criticism is that inflow periodic unsteadiness may trigger a diffuser instability, see e.g., [3]. Another criticism is that the theory assumes diffuser flow as inviscid. This implies that the dynamic character of diffuser rotating stall depends on the inviscid core flow and is independent of viscous effects, dominating in the boundary layer. It is important to acknowledge some observations where the periodic unsteadiness may be delaying stall, see e.g., [4] where low pressure turbine suction side boundary layer stall, and axial compressor stall were investigated.
For the actual evolution of the dynamic character and thus prediction of the number of rotating stall cells and their propagation speed around the annulus, different research groups have considered a time-evolving calculation assuming 2D-unsteady, inviscid and incompressible flow, see e.g., [5,6,7]. The effect of the blades on the flow is fed at the diffuser inlet boundary and a Dirichlet condition with static pressure is considered on the outlet surface. One observation is that the evolution of the rotating stall cells depends on the specified outlet boundary condition; see e.g., [8]. Specification of a Dirichlet boundary condition at the outlet was seen to result in strong gradients of the velocity field. This is in contrast to the linear stability analysis where the Neumann boundary condition is used, which assumes zero pressure gradients in both the radial and tangential directions. In scenarios with reversing flow at the outlet boundary special consideration is needed for an appropriate specification of the backflow direction. This could be extrapolated or it could be treated as normal to the face elements on the outlet surface. In the ideal case, the outlet boundary would be extended further downstream to prevent recirculation. In the linear stability analysis work of [1], extrapolating the backflow direction was seen to be in agreement with experimental results of [9] for a wide vaneless diffuser with uniform inflow. The rotating stall feature was illustrated as two high static pressure cells at low velocity with reversing flow and two lower pressure cells at higher outgoing velocity. This was later confirmed in the 2D linear stability analysis by [5]. These cells were reported to be rotating at sub-synchronous rotational speed. In a real diffuser design the channel width is narrow where viscous effects in the boundary layer cannot be neglected. This is due to the growth of the displacement thickness that may be different on the wall boundaries, depending on the quality of the flow delivered by the impeller, see e.g., [10]. Additionally, the rotation of the impeller with high blade tip speed results in a high Reynolds number flow that enters the diffuser. Therefore the flow is thus turbulent with large fluctuations in the flow quantities, and subsequently with high shear rate, which may become unstable. Typically, the fluctuations are amplified at off-design conditions, as documented by [11]. These effects are not taken into account with the 2D inviscid approach. In the work by [12] the Navier–Stokes equations were averaged with the Reynolds decomposition resulting in the need for turbulence closure modeling. They employed a two equation isotropic turbulence model. The numerical result was found in good agreement with observed pressure fluctuations and number of rotating stall cells as compared with particle imaging velocimety (PIV), see [13]. However, it was concluded that the numerical result using isotropic turbulence modeling underestimates the strength of the rotating stall vortex as compared with PIV measurements. This may be explained by excessive turbulent diffusion, a known issue of linear turbulence models in presence of adverse pressure gradients, and issues associated with predicting accurately features such as boundary layer separation. Many different closure models exist for calculating turbulent flow. In the case of the diffuser flow there is evidently some argument for including curvature effects in the modeling approach, since curvature requires discerning turbulence anisotropy. Consequently, it is clear that the turbulence model should take anisotropic effects into consideration on top of considering an approach able to resolve the large, energy containing flow structures. Streamline-curvature is possible to include in the framework of eddy-viscosity turbulence modeling. For this, a correction factor is introduced to the turbulent production term in the turbulent kinetic energy equation, see e.g., [14]. This factor depends on the strain-rate as well as the rotation-rate tensor including several proposed modeling constants. Due to the large freedom in choosing the modeling constants a small variation to proposed default values may potentially yield significant differences in the predicted turbulence production. Nevertheless, eddy-viscosity modeling with curvature correction showed good predictive quality compared to experimental data for a vaneless diffuser application during rotating stall conditions, as shown in [15].
Previously it was demonstrated that the angular momentum instability associated with curvature plays a key role in provoking large effects on the turbulent boundary layer, see work by [16]. The boundary layer is affected by the radial pressure gradient, radial curvature near the inlet and tangential straining as the flow cross section area increases. Since the flow is swirling it is, therefore, necessary to include 3D effects and the destabilizing tangential curvature.
The present paper aims to capture the evolution of the rotating stall cells and their propagation around the annulus by employing the large eddy simulation (LES)-based approach to a 3D annular vaneless diffuser configuration with a turnaround outlet. The possibility of capturing anisotropic turbulence effect in the boundary layer will help elucidate the variation of velocities and Reynolds stresses in the swirling diffuser flow. Especially in key areas close to separation and recirculation in order to investigate the effect on fluctuation levels. For a statistical assessment the flow field is quantified by means of flow mode decomposition techniques. A sensitivity study to the outlet boundary conditions imposed will be carried out. The numerical methodology is described in the next section, which is followed by presentation of results. The relevance of the obtained results to the scientific and designer community is included in a concluding summary section.

2. Numerical Methodology

The narrow annular vaneless diffuser with a turnaround outlet, shown in Figure 1, contains a converging section close to the impeller blades. The overall extent before the turnaround in the radial direction is 1.5 R 2 . The curvature of the turnaround is approximately 0.2 R 2 and the outlet is extruded to a location at 0.82 R 2 .
The effect of the impeller blade passing is provided as meridional and tangential velocity profiles, respectively. The 3D time-dependent large-scale turbulent motion responsible for the fluid mixing in the vaneless diffuser is computed using the LES approach. The computational domain is discretized with a curvilinear grid which constitutes the elements of the finite volume methodology employed. Via Taylor series expansion an approximate representation is obtained for the transient, convective and diffusive terms respectively. The temporal term is discretized with a second order scheme, which uses the solution from the current time step and from the previous two time levels. A hybrid second order upwind/central differencing scheme is used for the convective term. The diffusive flux term is treated with a second-order expression, which involves the cell center value, the neighboring cell center value, and including the diffusion flux at the interior face. The introduced dissipative truncation error is known to mimic Smagorinsky-type SubGrid scale (SGS) modeling, see works by [17,18]. As the grid size is refined the filter size for SGS modeling subsequently decreases. Hence, the contribution from SGS may therefore become very small as compared to other entities in the momentum balance equation. One may therefore disregard from the SGS term by setting it to zero, which is attributed as implicit LES with no explicit SGS. The validity of an implicit LES approach can be verified by evaluating the energy decay slope for a point exhibiting isotropic turbulence in the inertial subrange, i.e., 5 / 3 for a velocity fluctuating component. However, such condition is not expected for a narrow channel flow under pressure gradient and strong flow curvature yielding significant effect on the boundary layer development and with associated separation. Therefore, code validation is assessed alternatively by means of a grid dependency study, which is presented further below. Nevertheless, the solver has been used previously for calculating flow scenarios associated with centrifugal compressors at low mass flow rates and validations against experimental data were carried out, see [19,20,21].
A constant-density turbulent flow is assumed with a Reynolds number R e b 1 × 10 5 and inlet reference Mach number M r e f = 0.2 . The reference flow angle at the inlet boundary is relative to the tangential. It is computed in a point between the blades (i.e., θ / χ = 0.5 ) and at the diffuser channel midpoint between the shroud and hub side (i.e., n = 0.5 ):
α I N = arctan ( U m / U θ )
At the inlet boundary U m and U θ have variable distributions in both axial and tangential directions, with purpose to mimic a real velocity distribution from an upstream impeller. Figure 2 shows the meridional and tangential velocity profile distributions scaled with U t i p for two different reference inlet flow angle cases α I N = 9.5 and α I N = 5.5 .
For each blade tip width segment circumferentially around the inlet diffuser plane, U θ = U t i p and U m = 0 . This has the effect of modeling the blade blockage effect, which is presented as tangential distributions for two blade passages, see last row in Figure 2 for orientation purposes. In addition, the velocity distribution is made to rotate at the impeller angular velocity f I M P . Comparing peak values in Figure 2 it can be seen that the amplitude of the flow distortion entering the domain is in the order of 0.1 U t i p at the distinct blade passing frequency f B P F . Additional stochastic fluctuation is not considered, due to the overwhelming effect of the periodic unsteady incoming wakes described in Figure 2. At the outlet boundary, the pressure is fixed when the flow is outgoing. In case of backflow the pressure and flow direction is extrapolated from the outlet boundary. At the inlet the velocity components are specified and pressure is allowed to float. Therefore, the pressure is fixed in the scenario with pure outflow at the pressure outlet cell faces.
For validity of the numerical result a grid dependency study is performed on three different grids named; coarse, medium and fine (see also Figure 3). A grid refinement factor of two is used between grids. The wall resolution for the fine grid considered is: Δ n + ¯ = 0.5 , Δ m + ¯ = 30 , and Δ θ + ¯ = 30 . The normalized wall distances are computed as the surface average of the hub and shroud walls, respectively, between meridional stations m = 0 and m = 1 . For a fair comparison the time-step size is adjusted between grids so that the time averaged convective Courant number is unity. Approximately 30 through flow times were simulated for each case. Five through flow times were needed to establish a suitable initial turbulent flow field. Therefore, data corresponding to 25 through flow times are used for statistical analysis of the flow field, which is adequate to capture a low frequency instability in the order of 0.02 f B P F . The integrated pressure rise Δ C p is computed as the difference in area averaged pressure coefficient at meridional stations m = 0 and m = 0.8 , respectively. This parameter is evaluated with respect to the grid resolution in the meridional direction Δ m + ¯ . Figure 3 indicates that the solution (i.e., Δ C p ) approach an asymptotic limit value as the grid is refined. The Richardson 2nd order extrapolation shows that the fine grid resolution is marginally different compared to a hypothetical infinite grid resolution. It is also evident that the relative error with respect to the infinite grid solution drops monotonically two orders of magnitude for each refinement order. This result correlates with the chosen second order scheme and an expected truncation error of O ( Δ m + 2 ¯ ) .

3. LES Data Analysis

Based on the grid dependency assessment, the solution obtained with the fine grid resolution is chosen for further analysis. Figure 4 shows distributions on the side view of the time-averaged pressure, velocity components, and Reynolds stress levels, respectively, for the flow angle α I N = 9.5 . The Reynolds stresses are scaled as a fluctuation intensity with respect to the mean local velocity. This velocity is obtained at the midway distance between the shroud and the hub and along the meridional coordinate. A characteristic pressure increase with a relatively steep gradient is seen until the top of the turnaround. After the top of the turnaround the cross section area decreases. This leads to a flow acceleration and thus to a static pressure decrease (see C p scalar distribution in Figure 4). At the turnaround bend a small pressure gradient is present in the wall-to-wall normal direction with slightly lower pressure on the hub side compared to the shroud side.
The time-averaged meridional and tangential velocity profiles, respectively, show a strong shear layer close to the inlet, which is due to the jet-like shape of the imposed inlet boundary profile. In the straight section the velocity components gradually approach the shape of a fully turbulent developed flow profile with high fluctuation levels close to the walls. To help see this more clearly the time-averaged meridional and tangential velocity profiles, respectively, are plotted in Figure 5 for different stations along the meridional direction, between m = 0 . 05 and m = 0.8 . At meridional m = 0.6 , just prior the turnaround, the profiles approach an approximate self-similar solution.
Beyond m = 0.6 both profiles are seen to tilt with reduced through flow close to the hub wall, an effect caused by the turnaround bend. Note that the time-averaged meridional velocity has a more extreme slope and indicates tendency of an inflection point between m = 0.7 and m = 0.8 . This leads to increased fluctuations (root mean square) of this component ( u m ). The steepness close to the shroud and hub at m = 0.6 is compatible with the law of the wall. Therefore, wall modeling with lower resolution near the wall may lead to acceptable predictions. However, in the vicinity of the developing inflection point of the meridional velocity component there is evidence that wall modeling would not be acceptable and compromise prediction of potential boundary layer separation.
The fluctuation components of velocity in the meridional, the tangential as well as the co-variant direction, respectively, are depicted in Figure 4. A more detailed distribution of the fluctuation intensities is also depicted in Figure 5 for different meridional stations. Close to the inlet, peak fluctuation intensity is located close to the wall. In the middle of the channel, the peak values reduce and approach self-similarity. One explanation for this is the shape of the imposed boundary condition with the blades passing by creating high stress levels in the wake flow. Further downstream however, the intensity of the fluctuation reduces, which indicates turbulence decay. It is noted that the Reynolds stress anisotropy component exhibits a sudden growth at m = 0.7 close to the hub wall (see Figure 4). However, this follows by a rapid recovery when the flow enters the area of favorable pressure gradient towards the outlet. In general the tangential stress level is dominant but the meridional and anisotropic stress components are not insignificant. This clarifies the importance of employing a numerical approach that can capture anisotropic features in diffuser flow.
Assessment of the level of curvature requires a determination of the flow angle with respect to the tangent. This is presented in Figure 6 at different meridional stations.
The wall is said to be in a collateral state if the flow angle α = arctan ( U m / U θ ) is independent of the wall normal direction n. In other words, there would be a linear relationship between tangential and meridional velocities. It appears that this may hold in the straight diffuser section. A collateral state does not seem to hold past the turnaround at m = 0.7 . The figure shows evidence of a change in the flow angle as moving downstream with two different flow angles at the hub and shroud walls, respectively. In other words a smaller flow angle holds on the hub side as compared to the shroud side, and presenting a moderate gradient in between. This could be due to the concerted action of the turnaround bending effect and viscous effect near the wall, affecting the swirling motion of the core flow differently in the near wall regions associated with the hub and the shroud, respectively. In addition to this, the flow angle is seen to increase gradually towards the turnaround, which therefore gives a stabilizing effect. The increasing flow angle can also be linked with a loss of angular momentum due to viscous effect. In the area of the turnaround the flow angle decreases close to the hub side, which is associated with reduced through flow. This clearly leads to a destabilizing effect of the flow.
The presented results so far have only exposed statistics of the flow variables. A further assessment is carried out analyzing the power spectral density (PSD) of a monitored flow quantity of interest. Figure 7a shows PSD for the tangential velocity at m = 0 . 1 in the mid of the diffuser channel, i.e., n = 0.5 . A combination of broadband and narrowband content is observed in the spectra. The tonality of the blade passing frequency including several higher harmonics is correctly captured with the numerical approach. In the low frequency end of the spectra for both inlet reference flow angles considered significant tonalities are revealed. For the larger flow angle case α I N = 9.5 , one notable peak is observed at 0.04 f B P F . For the lower flow angle case α I N = 5.5 two significant low frequency tonalities can be observed, one at 0.03 f B P F and the other at 0.07 f B P F . The low frequency tonalities (i.e., 0.04 f B P F for the large flow angle case and 0.03 f B P F and 0.07 f B P F for the smaller flow angle case), correspond to flow modes with a periodic wave-like character. For the lower flow angle case the intensity level in the low to middle frequency range is seen to amplify. Moreover, the peak of the dominant narrowband feature, i.e., 0.07 f B P F for α I N = 5.5 as compared to 0.04 f B P F for α I N = 9.5 , is shifted to a higher frequency relative to the shaft speed. It should be noted that the sample length corresponds to approximately 25 cycles of the low frequency tonality. Therefore, these features appear as broad peaks rather than sharp peaks in the point spectra. However, this could be improved with a more generous sampling length and hence simulation elapsed time. Although, if a harmonic signal is Gaussian shaped, the variation is not exactly sinusoidal.
Up to this point the numerical data elucidate two possible mechanisms for onset of the low frequency instability, i.e., one shear layer close to the inlet and a tendency of an inflection point in the turnaround section. To correlate these mechanisms with the rotating instability the point-to-point cross power spectral density (CPSD) of the tangential velocity fluctuation is presented in Figure 7b. The signals in the two different points show strong correlation at the rotating instability (i.e., 0.04 f B P F for α I N = 9.5 ) as well as at the blade passing frequency. The phase shift at the rotating instability frequency is approximately χ . This means that the signal closer to inlet is ahead approximately one blade-to-blade passage. For the smaller flow angle case α I N = 5.5 spatial coherency is also found in the low frequency range. There, the peak at 0.07 f B P F is seen to be more dominant over the peak at 0.03 f B P F . Since their frequency ratio is not a whole number integer they are not harmonics of each other.
For a statistical characterization of the flow mode associated with the interesting narrowband instability a Fourier surface spectra analysis is performed of the velocity fluctuations u r and u θ for every point distributed on a meridional surface located midway between the hub and the shroud, which follows the diffuser channel curvature. The number of cycles of the low frequency mode shape included in the surface spectra computation is limited to four and the result is presented in Figure 8 for the two different flow angles considered.
For a physical interpretation the intensities of the flow modes should be observed as a superposition on the mean flow ( U i ) where the subscript i is a direction index. The real part of the flow mode spectra for the larger flow angle reveals two areas with negative intensity and two areas with positive intensity respectively, and distributed on opposing sides of each other. A positive radial intensity level is interpreted as outgoing flow and a negative level is a local inversion of the radial velocity component. The imaginary part (not included in the figure) is a phase shifted 90 representation of the real part in the rotation direction of the shaft (clockwise in the figure). With the real and imaginary parts added together u i = { u i } cos ( 2 π f t ) { u i } sin ( 2 π f t ) the shape of the mode is seen to rotate in the same direction as the impeller rotation. In other words, the flow mode can be described as a rotating instability consisting of two counter-rotating flow cells relative to each other. This qualitative description is in agreement with results of other research groups, e.g., [5,9]. The same feature can be attributed also to the lower flow angle (second and third columns) but with a larger number of rotating cells distributed evenly around the circumferential. The tonality at 0.03 f B P F for α I N = 5.5 (see Figure 7a) shows three disturbances (second column in Figure 8), whereas the other peak at 0.07 f B P F shows five rotating cells (third column in Figure 8).
One notable difference with the larger flow angle α I N = 9.5 as compared to the lower flow angle α I N = 5.5 , can be seen in terms of their most dominant low frequency tonalities, respectively. That is 0.04 f B P F (first column in Figure 8) as compared to f = 0.07 f B P F (second column in Figure 8). For the coherent mode shape at 0.04 f B P F the intensity of the tangential fluctuating mode shape is more dominant compared to the radial fluctuating mode shape, which is rather weak in intensity. For peak 0.04 f B P F in the lower flow angle case, this is reversed, i.e., the radial fluctuation dominates over the tangential. The reason to this can be attributed to type and location of the outlet boundary condition for flow case α I N = 5.5 . For this case a fraction of the cell face elements on the outlet were subjected to reversed flow, which means that those face elements can no longer be considered as a real outlet boundary. Under such situation the numerical solution may not be unique. Therefore, for the lower flow angle case there is a fraction of the outlet that violates the well-posedness criteria, see [22]. Strictly speaking, a solution to a problem that is not well-posed does not make sense, and so the numerical result for the lower flow angle case cannot be fully trusted. The remedy to this is to consider a repositioning of the outlet where the flow is purely outgoing. However, a common result is that the intensity level of the rotating disturbance reaches a peak at the top of the turnaround and with a growth point location prior to the turnaround.
The propagation speed of the rotating stall cells can be analyzed using a space-time cross correlation. For this the tangential velocity fluctuating component u θ is considered for all grid points, equidistant arranged along the circumferential direction at meridional m = 0.6 , between the hub and shroud midway in the diffuser channel. Among these points one is chosen as a reference located at twelve o’clock. Subsequently, the space-time cross correlation is computed as:
R ( r i , Δ t ) = < u ( x i , t ) u ( x i + r i , t + Δ t ) > < u ( x i , t ) > 2 < u ( x i + r i , t + Δ t ) > 2
The cross-correlation is then normalized and the result for flow angle cases α = 9.5 and α = 5.5 are presented in Figure 9.
The fluctuating tangential velocity component utilized in the space-time correlation is obtained by reconstruction of the most energetic Fourier surface spectra flow modes. The most energetic modes are determined from Figure 7. Thus, for flow angle case α I N = 9.5 the contributions from tonalities at 0.04 f B P F and f B P F , respectively, are used. For the lower flow angle case α I N = 5.5 , the contributions from 0.03 f B P F , 0.07 f B P F and f B P F are used. A coherent pattern is observed with several inclined lines exhibiting strong correlation. This is directly linked with the propagation of the rotating cells. The slopes of the inclined lines are approximately constant and correspond to the local convection speed. For the large flow angle case, there are two disturbances along the circumference, which correlates with the two rotating stall cell structures in Figure 8, (i.e., 0.04 f B P F for α I N = 9.5 ). For the smaller flow angle case there are thus five disturbances, since there are five dominant rotating stall cell structures. For α I N = 9.5 the local convection speed is approximately 50% of the impeller speed. Now, the slope of the characteristic is different for α I N = 5.5 , which is approximately 30% of the impeller speed. From Figure 7b and flow angle case α I N = 9.5 the tonality 0.04 f B P F is approximately one order larger than the f B P F tonality. Similarly, for flow angle case α I N = 5.5 , the tonality at 0.07 f B P F is the most relevant. Therefore, only one set of inclined lines are appearing in each two-point cross-correlations (Figure 9).
Another interesting observation in Figure 8 is that the radial and tangential velocity fluctuation distributions elucidate streaky elongated flow features, which propagate with the rotating stall mode. The length scale of the streaks is in the order of the diffuser channel width. Close to the inlet these features are located in the viscous sublayer rather in the core flow, which correlate with the amplified Reynolds stress levels close to the wall as shown in Figure 5. They can be identified as features carrying turbulent energy, which is an important characteristic of wall-bounded turbulence. Further out in the radial direction these stream-wise streaks grow with increased interaction and mixing with the core flow, see Figure 10. This partially motivates why the Reynolds stress distribution in the tangential direction becomes fuller in the mid of the diffuser channel.
Qualitatively, there is a resemblance with the stream-wise streaks found in the turbulent boundary layer in curved channel flow. However, in this case there is the added complication of an interaction with a rotating stall mode feature and an adverse pressure gradient. It is known that secondary flow instability can appear in a turbulent boundary layer over a curved channel wall. An important factor is when the boundary layer thickness is comparable to the radius of curvature. In such scenario, a centrifugal action creates a pressure variation across the boundary layer. This is said to lead to formation of longitudinal vortices [23]. In Figure 10, for the selected meridional section, for the larger flow angle case α I N = 9.5 , two major rotating secondary flow structures are seen with counter clockwise rotation. These are located closer to the shroud wall. The other major vortex in between have a clockwise rotation and is located closer to the hub wall. Consequently, just ahead of the clockwise rotating vortex u n ¯ is positive and just behind it u n ¯ is negative. A similar sequence of updraft and downdraft of the wall normal velocity component between the hub and shroud walls can also be seen for the small flow angle case α I N = = 5.5 . The onset of the instability can be estimated with the Görtler number, which is the ratio of centrifugal to viscous forces in the boundary layer. It is defined as:
G = U θ ν θ R 1 / 2
where U-bulk flow velocity, θ -momentum thickness, ν -kinematic viscosity, and R-curvature radius of the wall. For the considered diffuser geometry this number exceeds 0.3, which is the critical limit for occurrence of the Görtler instability. This is one possible explanation for the coherent streaky features found superimposed on the rotating stall mode (see Figure 8).

4. Conclusions

The computation of the mean velocity components in a narrow vaneless diffuser at low flow angles suggests a possible connection between the tendency of boundary layer separation in an adverse pressure gradient and the tendency of flow circulation and propagation into a rotating stall consisting of counter-rotating flow structures. This is in agreement with published literature on the subject, e.g., [1,5]. For the diffuser geometry considered in this study a tendency of an inflection point was detected in the time-averaged meridional velocity profile close to the hub wall at the top of the turnaround location. This correlates with a sudden reduction in the near-wall stress level in the near-wall region, which is indicative of approaching boundary layer separation.
An emerging feature of the rotating stall instability is the formation of counter-rotating flow structures. Such features of the flow causes unsteady effects on the turbulence structures residing in the boundary layer. One important perturbation mechanism is linked with the pressure gradient developed and the impact of the streamwise curvature. This conclusion is associated with the peak Reynolds stress level moving closer to the core flow and the response in the near-wall hub region, where flow blockage occurs. Another interesting feature is the characteristic turbulent streaks found in the turbulent boundary layer, which are being convected with the rotating stall instability. Those are turbulent energy carrying features. They were found to form longitudinal vortices with length scale in the order of the diffuser channel width. These features share resemblance with the secondary flow instability in turbulent boundary layers subjected to curvature effect.
A strong shear-layer is located close to the inlet, which is due to the considered shape of the imposed boundary condition. This correlates with high Reynolds stress levels and is another possible mechanism for the onset of the rotating instability which is seen to grow further downstream in the diffuser. The point-to-point cross-correlation was computed in order to determine the phase differences between the dominant flow modes, i.e., the rotating instability mode and the blade passing frequency mode. This is understood to be a convolution over a finite period with similar notion to a Fourier transform with assumption of periodicity. Consequently, all notions of causality are lost and one can no longer refer to driver and response signals, but only signals which are correlated. Therefore, it is challenging to say which signal drives which. However, the numerical result shows a strong coherence between the shear-layer instability close to the inlet and the developing instability prior to the turnaround location.

Acknowledgments

This work was conducted at the Competence Center for Gas Exchange (CCGEx). The Swedish National Infrastructure for Computing (SNIC) via HPC2N and PDC, the Parallel Computing Center at KTH, are acknowledged for providing the computing resources.

Author Contributions

E.S. performed all pre-, solve, post-processing. The type of post-processing was agreed and reviewed with M.M. and V.M.. M.G. and E.B. originally conceived the simplified numerical test set-up. This is the result of intense numerical testing that was discussed directly with E.S.. M.G. and E.B. also decided the boundary conditions and directed the team to investigate flow properties already detected in General Electric proprietary industrial experiments. M.M. and V.M. have a long and consolidated experience in scale-resolved simulations. They both supervised the runs and the post-processing, and worked with E.S. in the writing of the paper to make sure it had the necessary scientific rigor and industrial practical relevance.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

C p pressure coefficient 2 γ M r e f 2 ( P P r e f 1 )
P r e f reference pressure at the inlet (Pa)
M r e f reference Mach number at the inlet (-)
r , θ , z cylindrical coordinates defined in Figure 1
m , n meridional and wall normal coordinates defined in Figure 1
U r , U θ , U m mean radial, tangential and meridional velocities (m/s)
U t i p blade tip velocity (m/s)
u r , u θ , u m fluctuating velocities in radial, tangential and meridional directions (m/s)
f I M P , f B P F impeller angular velocity and blade passing frequency (Hz)
Nnumber of rotating flow cells or number of grid points
bdiffuser channel width at the inlet radius R 2 (m)
R 2 diffuser inlet radius (m)
χ blade-to-blade stagger angle (radian)

References

  1. Jansen, W. Rotating stall in a radial vaneless diffuser. J. Basic Eng. 1964, 86, 750–758. [Google Scholar] [CrossRef]
  2. Abdelhamid, A. Analysis of Rotating Stall in Vaneless Diffusers of Centrifugal Compressors. In Proceedings of the ASME 1980 International Gas Turbine Conference and Products Show, New Orleans, LA, USA, 10–13 March 1980. [Google Scholar]
  3. Madhavan, S.; Wright, T. Rotating stall caused by pressure surface flow separation on centrifugal fan blades. J. Eng. Gas Turbines Power 1985, 107, 775–781. [Google Scholar] [CrossRef]
  4. Seifert, A.; Barabi, A.; Wygnanski, I. Delay of airfoil stall by periodic excitation. J. Aircr. 1996, 33, 691–698. [Google Scholar] [CrossRef]
  5. Tsujimoto, Y.; Yoshida, Y.; Mori, Y. Study of vaneless diffuser rotating stall based on two-dimensional inviscid flow analysis. J. Fluids Eng. 1996, 118, 123–127. [Google Scholar] [CrossRef]
  6. Ljevar, S.; De Lange, H.; Van Steenhoven, A. Two-dimensional rotating stall analysis in a wide vaneless diffuser. Int. J. Rotating Mach. 2006, 56420. [Google Scholar] [CrossRef]
  7. Kalinkevych, M.; Shcherbakov, O. Numerical modeling of the flow in a vaneless diffuser of centrifugal compressor stage. ISRN Mech. Eng. 2013, 602384. [Google Scholar] [CrossRef]
  8. Dawes, W. A simulation of the unsteady interaction of a centrifugal impeller with its vaned diffuser: Flow analysis. J. Turbomach. 1995, 117, 213–222. [Google Scholar] [CrossRef]
  9. Wachter, J.; Rieder, M. Influence of Design on the Onset and Behavior of Rotating Stall in a single Stage Centrifugal Compressor. VDI-Berichte 1985, 572, 591–605. [Google Scholar]
  10. Satish, K.V.V.N.K.; Guidotti, E.; Rubino, D.T.; Tapinassi, L.; Prasad, S. Accuracy of centrifugal compressor stages performance prediction by means of high fidelity CFD and validation using advanced aerodynamic probe. In Proceedings of the ASME Turbo Expo 2013: Turbine Technical Conference and Exposition, San Antonio, TX, USA, 3–7 June 2013. [Google Scholar]
  11. Gaetani, P.; Persico, G.; Mora, A.; Dossena, V.; Osnaghi, C. Impeller-vaned diffuser interaction in a centrifugal compressor at off design conditions. J. Turbomach. 2012, 134, 061034. [Google Scholar] [CrossRef]
  12. Dazin, A.; Coutier-Delgosha, O.; Dupont, P.; Coudert, S.; Caignaert, G.; Bois, G. Rotating instability in the vaneless diffuser of a radial flow pump. J. Therm. Sci. 2008, 17, 368–374. [Google Scholar] [CrossRef]
  13. Dazin, A.; Cavazzini, G.; Pavesi, G.; Dupont, P.; Coudert, S.; Ardizzon, G.; Caignaert, G.; Bois, G. High-speed stereoscopic PIV study of rotating instabilities in a radial vaneless diffuser. Exp. Fluids 2011, 51, 83–93. [Google Scholar] [CrossRef]
  14. Shur, M.L.; Strelets, M.K.; Travin, A.K.; Spalart, P.R. Turbulence modeling in rotating and curved channels: assessing the Spalart-Shur correction. AIAA J. 2000, 38, 784–792. [Google Scholar] [CrossRef]
  15. Marconcini, M.; Bianchini, A.; Checcucci, M.; Ferrara, G.; Arnone, A.; Ferrari, L.; Biliotti, D.; Rubino, D.T. A Three-Dimensional Time-Accurate Computational Fluid Dynamics Simulation of the Flow Field Inside a Vaneless Diffuser During Rotating Stall Conditions. J. Turbomach. 2017, 139, 021001. [Google Scholar] [CrossRef]
  16. Clausen, P.; Koh, S.; Wood, D. Measurements of a swirling turbulent boundary layer developing in a conical diffuser. Exp. Therm. Fluid Sci. 1993, 6, 39–48. [Google Scholar] [CrossRef]
  17. Margolin, L.G.; Rider, W.J. A rationale for implicit turbulence modelling. Int. J. Numer. Methods Fluids 2002, 39, 821–841. [Google Scholar] [CrossRef]
  18. Fureby, C.; Grinstein, F.F. Large eddy simulation of high-Reynolds-number free and wall-bounded flows. J. Comput. Phys. 2002, 181, 68–97. [Google Scholar] [CrossRef]
  19. Semlitsch, B.; Mihaescu, M. Flow phenomena leading to surge in a centrifugal compressor. Energy 2016, 103, 572–587. [Google Scholar] [CrossRef]
  20. Sundström, E.; Semlitsch, B.; Mihaescu, M. Assessment of the 3D Flow in a Centrifugal compressor using Steady-State and Unsteady Flow Solvers. SAE Tech. Papers 2014. [Google Scholar] [CrossRef]
  21. Sundström, E.; Semlitsch, B.; Mihaescu, M. Generation Mechanisms of Rotating Stall and Surge in Centrifugal Compressors. Flow Turbulence Combust. 2017. [Google Scholar] [CrossRef]
  22. Hadamard, J. Sur les problèmes aux dérivées partielles et leur signification physique. Princeton Univ. Bull. 1902, 13, 28. [Google Scholar]
  23. Saric, W.S. Görtler vortices. Annu. Rev. Fluid Mech. 1994, 26, 379–409. [Google Scholar] [CrossRef]
Figure 1. (a) An isometric view of the vaneless diffuser geometry with a cut-out to depict the side view profile; (b) side section view depicting inlet and outlet stations.
Figure 1. (a) An isometric view of the vaneless diffuser geometry with a cut-out to depict the side view profile; (b) side section view depicting inlet and outlet stations.
Ijtpp 02 00019 g001
Figure 2. Scaled meridional and tangential velocity profiles applied at the diffuser inlet boundary condition for the two different flow angle cases considered. The graphs on the top row of the figure are for a station between the blades (i.e., θ / χ = 0.5 ). The graphs on the bottow row of the figure are for the midpoint between the hub and the shroud (i.e., n = 0.5 ). The small sketch to the right illustrates that the flow angle α is defined with respect to the tangent.
Figure 2. Scaled meridional and tangential velocity profiles applied at the diffuser inlet boundary condition for the two different flow angle cases considered. The graphs on the top row of the figure are for a station between the blades (i.e., θ / χ = 0.5 ). The graphs on the bottow row of the figure are for the midpoint between the hub and the shroud (i.e., n = 0.5 ). The small sketch to the right illustrates that the flow angle α is defined with respect to the tangent.
Ijtpp 02 00019 g002
Figure 3. Grid dependency assessment based the on Δ C p between m = 0 and m = 0.8 for three different spatial resolutions. The Richardson second order extrapolation indicates the hypothetical infinite grid solution. The relative error is seen to reduce two orders magnitude with respect to the grid size Δ m + ¯ .
Figure 3. Grid dependency assessment based the on Δ C p between m = 0 and m = 0.8 for three different spatial resolutions. The Richardson second order extrapolation indicates the hypothetical infinite grid solution. The relative error is seen to reduce two orders magnitude with respect to the grid size Δ m + ¯ .
Ijtpp 02 00019 g003
Figure 4. Time-averaged scalar contour distributions for α I N = 9.5 , colored on the side view. Flow quantities of interest are indicated in the text annotations.
Figure 4. Time-averaged scalar contour distributions for α I N = 9.5 , colored on the side view. Flow quantities of interest are indicated in the text annotations.
Ijtpp 02 00019 g004
Figure 5. Profiles of the time-averaged meridional and tangential velocity components for flow angle α I N = 9.5 for different locations along the meridional indicated in the figure. Profiles of mean meridional and tangential velocity fluctuation intensity are shown in the bottom row. Profiles are scaled with the mean local velocity.
Figure 5. Profiles of the time-averaged meridional and tangential velocity components for flow angle α I N = 9.5 for different locations along the meridional indicated in the figure. Profiles of mean meridional and tangential velocity fluctuation intensity are shown in the bottom row. Profiles are scaled with the mean local velocity.
Ijtpp 02 00019 g005
Figure 6. Profiles of the time-averaged flow angle relative the tangent for reference inlet flow angle α I N = 9.5 for different locations along the meridional indicated in the figure.
Figure 6. Profiles of the time-averaged flow angle relative the tangent for reference inlet flow angle α I N = 9.5 for different locations along the meridional indicated in the figure.
Ijtpp 02 00019 g006
Figure 7. (a) Power spectral density of the tangential velocity component in a probe point located at m = 0 . 1 in the mid of the diffuser channel n = 0.5 . Data is compared for the two different flow angles considered. (b) Cross power spectral density of the tangential velocity component between probe points located at m = 0.1 and m = 0.7 in the mid of the diffuser channel n = 0.5 . Data is for the two flow angle case considered.
Figure 7. (a) Power spectral density of the tangential velocity component in a probe point located at m = 0 . 1 in the mid of the diffuser channel n = 0.5 . Data is compared for the two different flow angles considered. (b) Cross power spectral density of the tangential velocity component between probe points located at m = 0.1 and m = 0.7 in the mid of the diffuser channel n = 0.5 . Data is for the two flow angle case considered.
Ijtpp 02 00019 g007
Figure 8. Fourier surface spectra flow mode decomposition. Narrowband radial (top row) and tangential (bottom row) velocity fluctuation intensity colored on a midway-section between hub and shroud which follows the diffuser curvature. In detail, this correspond to n = 0.5 , m [ 0 , 0.8 ] , θ [ 0 , 2 π ] . In the figure this is presented in a frontal view, which means that points for meridional positions at the top of the domain and beyond towards the outlet are not seen. Rotating stall modes at different flow angles, α I N = 9.5 (first column) and α I N = 5.5 (second and third columns). Last column shows the blade passing frequency flow mode for α I N = 9.5 .
Figure 8. Fourier surface spectra flow mode decomposition. Narrowband radial (top row) and tangential (bottom row) velocity fluctuation intensity colored on a midway-section between hub and shroud which follows the diffuser curvature. In detail, this correspond to n = 0.5 , m [ 0 , 0.8 ] , θ [ 0 , 2 π ] . In the figure this is presented in a frontal view, which means that points for meridional positions at the top of the domain and beyond towards the outlet are not seen. Rotating stall modes at different flow angles, α I N = 9.5 (first column) and α I N = 5.5 (second and third columns). Last column shows the blade passing frequency flow mode for α I N = 9.5 .
Ijtpp 02 00019 g008
Figure 9. Space-time cross correlation of the tangential velocity fluctuating component along the azimutal direction at m = 0.6 and in the middle of the diffuser section at n = 0.5 . (a) Flow angle α I N = 9.5 , tangential velocity fluctuation reconstructed from Fourier surface spectra modes at 0.04 f B P F and f B P F , respectively; (b) flow angle α I N = 5.5 , tangential velocity fluctuation reconstructed from modes at 0.03 f B P F , 0.07 f B P F and f B P F . The correlation is normalized, and the scale range is from −1 to 1.
Figure 9. Space-time cross correlation of the tangential velocity fluctuating component along the azimutal direction at m = 0.6 and in the middle of the diffuser section at n = 0.5 . (a) Flow angle α I N = 9.5 , tangential velocity fluctuation reconstructed from Fourier surface spectra modes at 0.04 f B P F and f B P F , respectively; (b) flow angle α I N = 5.5 , tangential velocity fluctuation reconstructed from modes at 0.03 f B P F , 0.07 f B P F and f B P F . The correlation is normalized, and the scale range is from −1 to 1.
Ijtpp 02 00019 g009
Figure 10. Fourier surface spectra decomposition of the wall normal velocity fluctuation colored on the side view plane, between meridional stations m = 0.3 and m = 0.4 . Constrained streamlines are overlayed for the corresponding to the flow modes α I N = 9.5 and α I N = 5.5 , respectively.
Figure 10. Fourier surface spectra decomposition of the wall normal velocity fluctuation colored on the side view plane, between meridional stations m = 0.3 and m = 0.4 . Constrained streamlines are overlayed for the corresponding to the flow modes α I N = 9.5 and α I N = 5.5 , respectively.
Ijtpp 02 00019 g010

Share and Cite

MDPI and ACS Style

Sundström, E.; Mihăescu, M.; Giachi, M.; Belardini, E.; Michelassi, V. Analysis of Vaneless Diffuser Stall Instability in a Centrifugal Compressor. Int. J. Turbomach. Propuls. Power 2017, 2, 19. https://doi.org/10.3390/ijtpp2040019

AMA Style

Sundström E, Mihăescu M, Giachi M, Belardini E, Michelassi V. Analysis of Vaneless Diffuser Stall Instability in a Centrifugal Compressor. International Journal of Turbomachinery, Propulsion and Power. 2017; 2(4):19. https://doi.org/10.3390/ijtpp2040019

Chicago/Turabian Style

Sundström, Elias, Mihai Mihăescu, Marco Giachi, Elisabetta Belardini, and Vittorio Michelassi. 2017. "Analysis of Vaneless Diffuser Stall Instability in a Centrifugal Compressor" International Journal of Turbomachinery, Propulsion and Power 2, no. 4: 19. https://doi.org/10.3390/ijtpp2040019

APA Style

Sundström, E., Mihăescu, M., Giachi, M., Belardini, E., & Michelassi, V. (2017). Analysis of Vaneless Diffuser Stall Instability in a Centrifugal Compressor. International Journal of Turbomachinery, Propulsion and Power, 2(4), 19. https://doi.org/10.3390/ijtpp2040019

Article Metrics

Back to TopTop