Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

55, 2, pp. 635-647, Warsaw 2017
DOI: 10.15632/jtam-pl.55.2.635

A NUMERICAL APPROACH TO PREDICT THE ROTATING STALL IN
THE VANELESS DIFFUSER OF A CENTRIFUGAL COMPRESSOR USING

THE EIGENVALUE METHOD

Chenxing Hu, Pengyin Liu, Xiaocheng Zhu

Co-Innovation Center for Advanced Aero-Engine, Shanghai Jiao Tong University, Shanghai, China

e-mail: ryanhu@sjtu.edu.cn; pengyinliu@163.com; zhxc@sjtu.edu.cn

Hua Chen

National Laboratory of Engine Turbocharging Technology, Tianjin, China

e-mail: huachen204887@163.com

Zhaohui Du

Co-Innovation Center for Advanced Aero-Engine, Shanghai Jiao Tong University, Shanghai, China

e-mail: zhdu@sjtu.edu.cn

Atwo-dimensional incompressibleflowmodel is presentedto study theoccurrenceof rotating
stall in vaneless diffusers of centrifugal compressors.The diffuser considered has two parallel
walls, and theundisturbedflow is assumedtobe circumferentiallyuniform, isentropic, and to
havenoaxial velocity.The linearized2DEuler equations for an incompressible flow in afixed
frame of the coordinate system are considered. After discretization by a spectral collocation
method based on Chebyshev-Gauss-Lobatto points, the generalized eigenvalue problem is
solved through the QZ algorithm. The compressor stability is judged by the imaginary part
of the eigenvalue obtained. Based on the 2D stability analysis, the influence of inflow angle,
radius ratio and wave number are studied. The results from the present stability analysis
are compared with some experimental measurement and Shen’s model. It is showed that
diffuser instability increases rapidly and the stall rotational speed decreases quicklywith an
increase in the diffuser radius ratio. The largest critical inflow angle can be obtained when
the wave number is around 3∼ 5 for the radius ratio between 1.5 to 2.2. It is also verified
that the stability model proposed in this paper agrees well with experimental data and has
the capability to predict the onset of rotating stall, especially for wide diffusers.

Keywords: instability, vaneless diffuser, eigenvalue problem, spectral method

1. Introduction

Rotating stall in radial vaneless diffuser is oneof themost commonflow instabilities in centrifugal
compressors and it can significantly influence the performance of the compressors. The nature of
flow instability, especially rotating stall associated with vaneless diffusers have been extensively
investigated bynumerous researchers (Day, 2016; Everitt andSpakovszky, 2013; Spakovszky and
Roduner, 2009; Ubben and Niehuis, 2015).

For decades, efforts have beenmade by researchers to explain themechanismand predict the
occurrence of rotating stall within the vaneless diffuser. Quantities of theoretical methodologies
based on different assumptions and simplifications on the base flow have been proposed du-
ring that time. Generally speaking, one of the basic distinctions between the theoretical models
on vaneless diffusers is whether the influence of the boundary layer is taken into considera-
tion. The first type of approaches which was adopted by Jansen (1964), Senoo and Kinoshita
(1977) and Frigne and Braembussche (1985) is that the three-dimensional wall boundary layer



636 C. Hu et al.

is supposed to be responsible for the occurrence of rotating stall in a narrow vaneless diffuser.
According to the experimental and theoretical study of Jansen (1964), the flow was assumed
to be symmetric with respect to the diffuser depth and the local inward radial velocity com-
ponent was treated as the inception of rotating stall. While in Senoo’s vaneless diffuser model,
the flow was no longer assumed to be symmetric and a non-uniform distribution of inlet velo-
city along the axial direction was taken into account. In his calculations, the critical velocity
angle was defined for which reverse flow started in the vaneless diffuser (Senoo and Kinoshita,
1978).

On the other hand, an other theory for the occurrence of rotating stall in a vaneless diffu-
ser, which was employed by Abdelhamid (1980) and Moore (1989), indicates that the stall is
associated with two-dimensional core flow instability in vaneless diffusers which usually have
the width radius ratio above 0.1. The two-dimensional numerical model developed by Moore is
based on calculation of 2D incompressible Euler’s equations. Neutrally stable rotating distur-
bances with low speed were found and a dense set of resonant solutions in which large pressure
perturbations were taken as the criterion of rotating stall in Moore’s work. Chen et al. (2011)
extended the 2 dimensional model of Moore into a 3 dimensional model with consideration of
the distribution of inlet velocity along the axial direction. Sun et al. (2013, 2016) proposed sta-
bility models for axial and centrifugal compressors on the basis of the eigenvalue approach. The
comparisons with the results from experiments validated the effectiveness and accuracy of their
models.

In addition, the significant influence of diffuser geometry and flow parameters on the va-
neless diffuser performance and structure of the stall pattern have been also numerically and
experimentally investigated in the recent years. The experiments carried out by Abdelhamid
(1983) and Bianchini et al. (2013) confirmed that the local reverse flow did not necessari-
ly lead to stall. The dependency of flow instability on the diffuser width ratios and diame-
ter were also verified. The critical flow coefficient becomes larger with a decrease in diffu-
ser width according to experiments of Abidogun (2002). Besides those experimental resear-
ches, the diffuser stability was also numerically investigated by Everitt (2010) through con-
ducting isolated diffuser simulations, and the volute was found to have potential in delaying
the onset of diffuser instability. Although the CFD method has an advantage over theoreti-
cal methods by providing a direct and vivid flow field with a relatively high accuracy, there
is no a certain way for numerical simulation to capture the complicated disturbance with dif-
ferent frequency, amplitude and length scale due to unsteadiness and complexity of the flow
field.

According to Jansen (1964), unsteady inviscid motion of the fluid is analyzed with the
assumption that the disturbances can be expressed in terms of periodic waves. The equations
are then reduced and solutions are sought for the resulting eigenvalue problem. However, it is
mentioned in his paper that the prediction of the number of stall cells and the determination of
the constant in the velocity equation requires a subsequentstudy.

In the present paper, firstly a 2D vaneless diffuser theoretical model based on stability ana-
lysis of an incompressible base flow is established and the full eigenvalue spectrum is obtained.
Then the influence of collocation points and geometric parameters of the vaneless diffuser on the
stability is investigated. Finally, the comparison between the results obtained from the present
stability model and those from several experiments and models reported previously is perfor-
med. Compared to unsteady CFD simulations which are quite time and resources consuming
and Senoo’s stability model, which is unable to provide the stall number, the present analy-
sis is able to predict the occurrence and stall number cheaply and quickly. And compared to
Shen’smodel, the presentmodel is easier to be extended to the stability problem concerning the
sensitivity.



A numerical approach to predict the rotating stall... 637

2. Theoretical model

2.1. Numerical methodology

In the present analysis, the stability discussed here is assessed by linearizing 2D Euler’s
equations for a base flow and determined by an operator that describes the evolution of small
perturbations superposed on the base flow. The eigenvalues of this operator give the frequency
and time growth rate of the perturbations. The corresponding eigenfunctions yield the mode
shapes. The perturbationwith the largest growth rate is the one that dominates the stability in
the long term, and the positive growth rate means the occurrence of instability.
The diffuser considered has two parallel walls, and the undisturbed flow is assumed to be

circumferentially uniform, isentropic, and to have no axial velocity. To further simplify the
calculation, the perturbations and velocity distribution of the base flow along the axial direction
are neglected. Theflow sketch is shown inFig. 1. The inflow circumferential angleα is defined at
the diffuser inlet as illustrated in Fig. 1b. Then, a 2D flowmodel for study of the flow stability
in the vaneless diffuser can be developed.

Fig. 1. Sketch of a vaneless diffuser: (a) base flow, (b) inflow angle

2.2. Implementation of the numerical model

2.2.1. Linearization of Euler’s equations

To begin with, viscosity and compressibility are neglected in our study. The flow field is
described by non-dimensionalized 2-D unsteady incompressible Euler’s equations as shown

vr
r
+
∂vr
∂r
+
1

r

∂vθ
∂θ
=0

∂vr
∂t
+vr
∂vr
∂r
+
vθ
r

∂vr
∂θ
−
v2θ
r
=−
∂p

∂r
∂vθ
∂t
+vr
∂vθ
∂r
+
vθ
r

∂vθ
∂θ
+
vθvr
r
=−
∂p

r∂θ

(2.1)

Here, the velocity is non-dimensionalized by the impeller tip velocity U. The length is non-
dimensionalized by the corresponding inlet radius r1 and the time t is non-dimensionalized
by r1/U. Thus, vr and vθ mean the non-dimensionalized radial and circumferential velocities. In
the stability analysis, theflowfield is assumed to consist of thebaseflowanda small disturbance.
Then, the velocity and pressure can be rewritten as follows

p= p+p′ vr = vr+v
′

r vθ = vθ+v
′

θ (2.2)



638 C. Hu et al.

wherep represents baseflowpressure and p′ are small disturbances of pressure.For thebaseflow,
Vr and Vθ represent the inlet radial and circumferential velocity of the base flow, respectively.
BeforeEqs. (2.2) are substituted intoEqs. (2.1), emphasis on the small perturbation theory is

made: 1) any disturbances of the 2nd and higher orders are neglected, 2) the base flow is treated
as steady, so ∂•/∂t=0, c) since the mean flow is circumferentially symmetric, ∂•/∂θ=0. And
Eqs. (2.1) are still appropriate for the base flow. Then the linearized non-dimensional Euler’s
equations can be derived as follows

v′r
r
+
∂v′r
∂r
+
1

r

v′θ
θ
=0

∂v′r
∂t
+
∂vr
∂rv′r

+
(
vθ
∂v′r
r∂θ
+vr
∂v′r
∂r

)
−
2vθv

′

θ

r
+
∂p′

∂r
=0

∂v′θ
∂t
+v′r

(∂vθ
∂r
+
vθ
r

)
+
(
vθ
∂v′θ
r∂θ
+vr
∂v′θ
∂r

)
+
vrv
′

θ

r
+
∂p′

r∂θ
=0

(2.3)

2.2.2. Establishment of the eigenvalue problem

In general, the separation between spatial and temporal coordinates allows making use of
Fourier modes in time when the variables are homogeneous. Then the perturbations can be
written in form of eiθ where the homogeneous variables (including time) are taken into account
and the inhomogeneous variables are taken as the amplitude function. The specific form of the
Fouriermodeassumed is that thedisturbances arenormalmodes in the circumferential direction,
the disturbances can be written as follows

p′ =P(r)ei(−ωt+mθ) v′r =A(r)e
i(−ωt+mθ) v′θ =W(r)e

i(−ωt+mθ) (2.4)

wherem is the azimuthal wavenumber and ω=ωr+iωi. Taking p
′ for instance as

p′ =P(r)eimθe(−iωr+ωi)t (2.5)

it can be seen that corresponding to eωit dominates the growth rate of disturbances with time.
Namely, if the imaginary part ωi > 0, the disturbances will grow exponentially with time, and
the mode will be unstable. If ωi < 0, the disturbances will decay with time and the mode will
be stable. And is the angular frequency.
The rotational speed of the disturbances in the circumferential direction non-dimensionlized

by the impeller rotational speed is

f =
−ωr
m

(2.6)

Substituting equations (2.4) into (2.3), the ordinary differential equations for A, W , P
which represent A(r), W(r) and P(r) with the subscript omitted, can be obtained as follows
(A′r = dA/dr etc.)

A′r +
A

r
+
im

r
W =0

− iωA+
∂vr
∂r
A+vrA

′

r+
imvθA

r
−
2vθ
r
W +P ′r =0

− iωW +A
(∂vθ
∂r
+
vθ
r

)
+vrW

′

r+
imvθW

r
+
vr
r
W +

im

r
P =0

(2.7)

Since there is no inlet disturbances coming from the outside of the system, homogeneous
ordinary differential equations (2.7) can also expressed as

Mφ=ωJφ (2.8)



A numerical approach to predict the rotating stall... 639

whereφ= [A,W,P]T, andM and J are the corresponding coefficient matrices

M=




1

r
+
∂

∂r

im

r
0

dvr
dr
+vr

∂

∂r
+
imvθ
r

−
2vθ
r

∂

∂r
dvθ
dr
+
vθ
r

vr
∂

∂r
+
imvθ
r
+
vr
r

im

r




J=



0 0 0
i 0 0
0 i 0


 (2.9)

Solving equation (2.8) leads to a generalized eigenvalue problem, and the eigenvalue is used
to determine the instability of the vaneless diffuser.

2.2.3. Establishment of the base flow

According to conservation of angular momentum and mass, the circumferential and radial
velocity can be expressed as follows

vr12πr1 = vr2πr vθ1r1 = vθr (2.10)

Therefore, the expressions of circumferential and radial velocity can be obtained as

vr =
Vr
r

vθ =
Vθ
r

(2.11)

where Vr and Vθ are the radial and circumferential velocity distributions at the diffuser inlet.

2.3. Discretization and calculation

2.3.1. Spatial discretization

Among the widely used spatial discretization methods, spectral methods have played a pro-
minent role in early instability analyses.At present they are particularly useful. In this study, the
spectralmethod is accomplishedwithChebyshev-Gauss-Lobatto points.The two-dimensional li-
near eigenvalue problems posed above can be solved by constructing an operator that discretizes
a n-th order linear ODE in the form of

Lij = an(ri)D
(n)
ij + . . .+a1(ri)D

(1)
ij +a0(ri) (2.12)

where D
(m)
ij , m = 1,2, . . . ,n, is the m-th derivative matrix corresponding to the collocation

points ri, and a0,a1, . . . ,an are evaluated at ri.

When using spectral collocation methods, the value of the functions at collocation points is
expressed as follows

φ(r)=
N∑

j=0

ϕj(r)φ̃(ri) (2.13)

whereφ(r) is the value of the function at the point r obtained by using interpolant polynomials
constructed for the variables in terms of their values at the collocation points, andϕ(r) is known
as the basis function.

The collocation points in the calculation domainΩ are chosen as

rΩ =cos
πj

N
j=0,1, . . . ,N (2.14)



640 C. Hu et al.

The extrema of the N-th order Chebyshev polynomial TN are defined in the interval
−1<rΩ < 1. Then the interpolant ϕj(rΩ) for the Chebyshev scheme is given as

ϕj(rΩ)=
1−r2Ωj
r−rΩi

T ′N(rΩ)

N2cj
(−1)j+1

c0 = cN =2 cj =1 0<j<N

(2.15)

The collocation derivative matrix D for the Gauss-Lobatto grid is denoted by

D
(1)
GL =(dij) 0¬ i j¬N (2.16)

where the elements d are defined by Canuto et al. (2010) as follows

d00 =−dNN =
2N2+1

6
djj =






rΩj
2(r2
Ωj
−1)

1¬ j¬N −1

ci
cj

(−1)i+j

rΩi−rΩj
i 6= j

(2.17)

In our case, the physical domain of interest r ∈ [1,Rf] is mapped onto the standard collo-
cation domain, based on the following equation

rΩ =−1+2
r−1

Rf−1
(2.18)

Thereby, linearized Euler’s equation (2.8) can be illustrated in the collocation domain as follows

MΩφ=ωJΩφ (2.19)

For the 2D incompressible problem, whenN collocation points are adopted to discretize the
domain of interest, MΩ is a matrix of size 3N × 3N arising from the three equations for the
problem.

2.3.2. Numerical calculation of the generalized eigenvalue problem

Onewidely used algorithm for solving generalized eigenvalue problems is theQZalgorithmof
Moler and Stewart (1973). This is a generalization of the QR algorithm for standard eigenvalue
problems.Comparedwith other algorithms such as theAnorldimethod (Meerbergen andRoose,
1997), inverse iteration (Peters andWilkinson, 1979) and the Jacobi-Davidsonmethod (Sleijpen
andVan derVorst, 1994), an important feature of theQZalgorithm is that it functions perfectly
well with the presence of an infinite eigenvalue due to singularity of thematrix jΩ. In the current
study, the calculation of the generalized eigenvalue problem proposed in equations (2.17) is
carried out by the QZ algorithm. Matlab is a software with high capability on computational
mathematics, and the function eig based on the QZ method in Matlab has been widely used
for eigenvalue problems such as for the calculation in the present model. For a matrix with
leading dimension of 200 in the presentmodel, the calculation with the function eig can be quite
time-saving.

What should be paidmuch attention to is that spurious eigenvalues derived from numerical
calculation rather than real physical problem also exist while using the QZ algorithm. Oneway
to distinguish the spurious eigenvalues from the real ones is to vary the number of collocation
points. The eigenvalues which remain the same are the real ones.



A numerical approach to predict the rotating stall... 641

3. Validation of the model

3.1. Influence of the collocation point number on the calculations

Firstly, the number of collocation points and its independence of the calculation is verified
through changing this number. The eigenvalue spectra with different numbers of collocation
pointsN at Rf =2, Vr =0.1, m=1 are shown in Fig. 2.
It can be seen that although the largest eigenvalue with the same imaginary part can be

obtained for all numbers of collocation points, the number larger than 50 is preferred to have a
proper, whole sketch of the eigenvalue spectrum. In the subsequent study, 50 collocation points
are adopted in the following calculations consideringbothnumerical accuracy andcomputational
efficiency. When the collocation number is more than 200 in the subsequent study, a 600×600
matrix is solved, and spurious eigenvalues which are at least two orders of magnitude higher
than the real ones come into being.

Fig. 2. Eigenvalue spectrum for different numbers of collocation points

3.2. Influence of diffuser inlet velocity of the base flow on the calculations

Generally, the inlet radial velocity distribution of the diffuser is highly influenced by the
impeller, and the axial uniform radial velocity can seldom be obtained in a practical diffuser
flow. However, the uniform distributed radial velocity for an inviscid flow is associated with the
average inflow angle. A brief connection between the average inflow angle and diffuser stability
can be revealed rapidly with the help of the present 2D stability model. Figure 3 shows the
influence of variability of the inlet velocity of the base flow Vr on the stability of the vaneless
diffuser with the same number N of collocation points, radius ratio Rf and wave number m.
Different inflow angles are also listed corresponding to different inlet velocities.



642 C. Hu et al.

It can be seen that with a decrease in the inlet radial velocity of the base flow, the flow in
the vaneless diffuser turns from stability to instability. In other words, with a decrease in the
inlet circumferential flow angle, the flow in the vaneless diffuser tends to instability, which well
agrees with Senoo and Kinoshita (1978). Correspondingly, the inflow circumferential velocity
angle when the rotating stall occurs is taken as the critical inflow angle αc.

Fig. 3. Eigenvalue spectrum for different inlet radial velocities of the base flow

3.3. Influence of geometric parameters on the instability

3.3.1. Influence of the radius ratio Rf

As shown inTable 1, six leading eigenvalues in a descending order of the imaginary part ofω
are presented.Comparedwith the cases forRf =1.5 and 1.67, the imaginary part of the leading
eigenvalue forRf =2 and 2.2 become positive. The leading eigenmode of the cases with larger
values ofRf turns to instability. A larger value of the imaginary part of the leading eigenvalue
indicates a higher time growth rate, whichmeans stronger instability of the flowfield.Therefore,
the flow tends to instability with an increase inRf.

3.3.2. Influence of the wave number m

A common feature of the centrifugal compression system is that it exhibits two kinds of
rotating stall. One is known as the impeller stall, having a relatively high rotational speed of
usually more than 50% of the impeller speed; and the other known as the diffuser stall, having
a slower rotational speed at about 10% of the impeller speed. The present stability model is
also applicable to confirmation of the stall number and rotational speed, therefore 3 cases with
different radius ratios Rf are calculated for 5 different wave numbersm.



A numerical approach to predict the rotating stall... 643

Table 1. Six leading eigenvalues for different values ofRf with Vr =0.131,N =50 andm=1

Rf =1.5 Rf =1.67 Rf =2 Rf =2.2

−0.10977−0.25458i −0.024169−0.05096i 0.03898+0.079355i 0.05377+0.104624i

−2.91568−0.64371i −2.22947−0.323566i −2.01654−0.236903i −0.29436−0.223245i

0.70504−0.72005i 0.332469−0.474960i −2.87547−0.395488i −1.7518−0.28062i

−4.28143−0.86662i −3.180583−0.49452i 0.327611−0.416233i −2.18808−0.32481i

2.04559−0.91810i 1.25444−0.59218i −3.73096−0.495389i −2.623871−0.33961i

−5.62630−1.01110i −4.12054−0.602551i 1.161860−0.514977i −0.08285−0.38274i

Whenm=0, the flow tends to be stable under any inlet flow condition. This is because that
the time growth rate of circumferential disturbances are themain concern in our study, and zero
wave number in the circumferential direction indicates that the circumferential disturbances are
neglected, resulting in the absolutely stable condition in the calculation.
As shown in Fig. 4a, the influence of the wave number m on the critical stall angle is

illustrated. It is clearly shown that the critical angle decreases with a decrease inRf. According
toAbidogum (2002), themost unstablemode (or stall) is often obtained for thewave numberm
of 3. It can be seen that the largest critical angle is obtainedwhenm is around 3with the radius
ratio of 2.2. This is consistent with the result of Abidogum. Thewave number corresponding to
the most unstable mode differs for various radius ratios.
In Fig. 4b, the rotational speed of the disturbance corresponding to the onset of instability f

vs. the wave number m is illustrated. It can be seen that the rotational speed of disturbance
linearly increases with an increase in thewave numbermwhile thewave speed decreases quickly
with the diffuser radius ratio, which is consistent with the results of Moore (1989) and Chen et
al. (2011).

Fig. 4. Influence of the wave numberm on the critical angle and rotational speed: (a) critical angle vs.
wave number, (b) rotational speed vs. wave number

3.4. Some eigenfunctions

In the theoryof linear stability, themaximumresponsetovaryiable initial conditionsarewhat
we concern themost. In the present analysis, the amplification of disturbances is characterized in
terms of the least stable mode of thematrixM. The eigenvalue spectrumwhen the radius ratio
Rf is 2, wave numberm is 1 and the inlet radial velocity is 0.04 are shown in Fig. 5. Among the
eigenvalues obtained, the leading three least stable eigenvalue spectra are selected and labeled as
ω1,ω2,ω3. The eigenfunctions for pressure andvelocity perturbations corresponding to the three
different eigenvalues are illustrated in Fig. 6, respectively. The eigenfunctions corresponding to
the three eigenvalues areproved tobe sortof non-orthogonal.Thenon-normality of the linearized
Euler equations has influence on the growth rate of the perturbations in a short timescale.



644 C. Hu et al.

Fig. 5. Eigenvalue spectrum forRf =2,m=1, V
r
=0.04

Fig. 6. Eigenfunctions corresponding to three eigenvalues for pressure and velocity perturbations:
(a) – ω1, (b) – ω2, (c) – ω3

4. Application of the present stability model

4.1. Comparison with experimental results in (Kinoshita and Senoo, 1985)

In order to verify the results of the present analysis, the experiment results byKinoshita and
Senoo (1985) are comparedwith the predictions from the 2D stability model. In the experiment
by Kinoshita and Senoo (1985), a back sweep angle of β = 24.6◦ was set at the impeller exit,



A numerical approach to predict the rotating stall... 645

and two different values of the radial ratio Rf were tested. The experiment was undertaken at
a low speed, and the inlet velocity distribution was not reported. The present 2D incompres-
sible stability model has been applied to the experimental set up with a uniform inlet velocity
distribution.
Since the rotating stall number observed in (Kinoshita andSenoo, 1985) is 3, the comparison

of the critical flow angle in (Kinoshita and Senoo, 1985), Shen’s in (Chen et al., 2011) and the
presentmethod for the wave numberm=3 is given in Table 2. It can be shown that the result
of the present stability analysis is much closer to the experiential data by Kinoshita and Senoo
(1985) than those derived from Shen’s 3D model.

Table 2. Comparison of the critical flow angle from experiment (Kinoshita and Senoo, 1985),
Shen’s result and current analysis

Rf β [◦] m
α [◦]

experiment Shen present model

1.4 24.6 3 5.2 7.925 6.657

1.67 24.6 3 5.5 12.82 11.631

At a smaller Rf, the agreement between the current method and the experiment is good.
At the same time, this study suggests that the flow becomes significantly less stable for a larger
radial ratio, whereas the experiment shows only a small deterioration of the stability.

4.2. Comparison with the experimental results in (Abidogun, 2002)

In (Abidogun, 2002) the experimental rig, as shown in Fig. 7, consists of an impeller with
radial blades and a vaneless diffuser. Three different values ofRf were tested. According to the
results given in (Abidogun, 2002), the stall wave numberwere always between 2-4.Whatmakes
rig in (Abidogun, 2002) different from that in (Kinoshita and Senoo, 1985) is the larger axial
width bz (> 0.1), which ismore suitable for themodels based on inviscid core flow theories such
as Shen’s and the present analysis.

Fig. 7. The experimental rig for a vaneless diffuser in (Abidogun, 2002)

It is illustrated in Table 3 that the inlet critical flow velocity at the onset of stall derived
by Abidogun (2002), Shen’s model and the present stability analysis shows the same trend of
increase with the diffuser radial ratio Rf. For a smaller radius ratio of the diffuser Rf, the
prediction results of the presentmodel and Shen’smodel show a higher accuracy than those for
the radius ratio 2.0. Compared with the results in (Kinoshita and Senoo, 1985), the predicted
critical inlet velocity for different radius ratios obtained in this test are all much closer to the
experimental results. Onemain reason for this lies in that the present analysis and Shen’smodel
are both based on the core flow assumption, which is more suitable for vaneless diffusers with
large width.



646 C. Hu et al.

Table 3. Comparison of the inlet critical flow in (Abidogun, 2002), Shen’s and the present
stability analysis

Rf m
Vrc

Abidogun Shen present model

1.5 3 0.16 0.1525 0.152

1.7 3 0.216 0.2150 0.215

2.0 3 0.255 0.2925 0.294

5. Conclusions

In this research, a twodimensional stabilitymodel of a vaneless diffuser in a centrifugal compres-
sor based on stability analysis is proposed. The independence of the number collocation points
and the influence of geometric parameters on the stability of the vaneless diffuser are verified.
The capability of the present model to predict the onset of stall with different radius ratios is
validated against several experimental data. The present model is effective in stall prediction
for a wide vaneless diffuser. Compared toMoore’s and Shen’s stability models, the advantage of
the present model lies in direct calculation and providing both the complex eigenspectrum and
rotational speed instead of multiple trials and iterations. Another remarking superiority is that
the present model investigates the instability induced by the inviscid main flow, and provides
the fundamental research for further study of the unsteady interaction of the inviscidmain flow
and the boundary layers.

In our investigation, the significant effects of the diffuser radius ratio on diffuser stability are
confirmed. The stability deteriorates rapidly with an increase in the radius ratio. The largest
critical inflow angle is obtained when the wave number m is around 3-5 for the radius ratio
between 1.5 to 2.2. Through model assessment, this model has the capability of predicting the
onset of stall in vaneless diffusers and can be applied in the cases without considering the axial
distribution of inlet velocity. However, in the cases where the axial distribution of inlet velocity
plays significant role, the application of a 3D stability model is required.

References

1. AbdelhamidA.N., 1980,Analysis of rotating stall in vaneless diffusers of centrifugal compressors,
Canadian Aeronautics and Space Journal, 26, 118-128

2. Abdelhamid A.N., 1983, Effects of vaneless diffuser geometry on flow instability in centrifugal
compression systems,Canadian Aeronautics and Space Journal, 29, 259-266

3. Abidogun K.B., 2002, Effects of vaneless diffuser geometries on rotating stall, Journal of Propul-
sion and Power, 22, 3, 542-549

4. BianchiniA.,BiliottiD.,BelardiniE.,GiachiM.,TapinassiL.,VanniniG.,FerrariL.,
FerraraG., 2013,Asystematicapproachto estimate the impactof theaerodynamic force induced
by rotating stall in a vaneless diffuser on the rotordynamic behavior of centrifugal compressor,
Journal of Engineering for Gas Turbines and Power, 135, 11

5. CanutoB.C., HussainiM.Y.,Quarteroni A., ZangT.A., 2010,SpectralMethods: Evolution
to Complex Geometries and Applications to Fluid Dynamics, Springer

6. Chen H., Shen F., Zhu X.C., Du Z.H., 2011, A 3D compressible flowmodel for weak rotating
waves invanelessdiffusers.Part I:Themodel andMachnumber effects,Journal ofTurbomachinery,
134, 4, 041010

7. Day I.J., 2016, Stall, surge, and 75 years of research, Journal of Turbomachinery, 138



A numerical approach to predict the rotating stall... 647

8. Everitt J.N., 2010, Investigation of stall inception in centrifugal compressors using isolated dif-
fuser simulations,M.S. Thesis, MIT, Cambridge,MA

9. Everitt J., Spakovszky Z., 2013, An investigation of stall inception in centrifugal compressor
vaned diffuser, Journal of Turbomachinery, 135

10. Frigne P., Braembussche R.V.D., 1985, A theoretical model for rotating stall in the vaneless
diffuser of a centrifugal compressor, Journal of Engineer for Gas Turbines and Power, 107, 2.
507-513

11. JansenW., 1964,Rotating stall in a radial vaneless diffuser,ASMEJournal of Basic Engineering,
86, 4, 750-758

12. Kinoshita Y., Senoo Y., 1985, Rotating stall induced in vaneless diffusers of very low specific
speed centrifugal blowers, Journal of Engineering for Gas Turbines and Power, 107, 514-521

13. Meerbergen K., Roose D., 1997, The restarted Arnoldi method applied to iterative linear
system solvers for the computation of rightmost eigenvalues, SIAM Journal on Matrix Analysis
and Applications, 18, 1, 1-20

14. Moler C.B., Stewart G.W., 1973, An algorithm for generalized matrix eigenvalue problems,
SIAM Journal on Numerical Analysis, 10, 2, 241-256

15. Moore F.K., 1989,Weak Rotating flow disturbances in a centrifugal compressor with a vaneless
diffuser, Journal of Turbomachinery, 111, 442-449

16. Peters G., Wilkinson J.H., 1979, Inverse iteration, Ill-conditioned equations and Newton’s
method, SIAM Review, 21, 3, 339-360

17. Senoo Y., Kinoshita Y., 1977, Influence of inlet flow condition and geometries of a centrifugal
vaneless diffuser on critical flow angle for reverse flow, Journal of Fluids Engineering, 99, 1, 98-103

18. SenooY.,KinoshitaY., 1978, Limits of rotating stall and stall in vaneless diffuser of centrifugal
compressors,ASME Paper 78-GT-19

19. Sleijpen G.L.G., Van der Vorst H.A., 1994, A Jacobi-Davidson iteration method for linear
eigenvalue problems, SIAM Review, 17, 2, 267-293

20. SpakovszkyZ., RodunerC., 2009, Spike andmodal stall inception in an advanced turbocharger
centrifugal compressor, Journal of Turbomachinery, 131

21. SunX.F.,LiuX.H., SunD.K., SunX., 2013,Atheoreticalmodel for flow instability inception in
transonic fan/compressors,ASMETurbo Expo 2013: Turbine Technical Conference and Exposition

22. Sun X.F., Ma Y.F., Liu X.H., Sun D.K., 2016, Flow stability model of centrifugal compressor
based on eigenvalue approach,AIAA Journal

23. Ubben S., Niehuis R., 2015, Experimental investigation of the diffuser vane clearance effect in
a centrifugal compressor stage with adjustable diffuser geometry. Part I: Compressor performance
analysis, Journal of Turbomachinery, 137

Manuscript received July 5, 2016; accepted for print December 17, 2016