Stability and structural sensitivity analysis of the turbulent flow in the narrow vaneless diffuser with mean flow method

Accepted Manuscript

Stability and structural sensitivity analysis of the turbulent flow in the narrow vaneless diffuser with mean flow method Chenxing Hu , Xiaojian Yang , Xiaocheng Zhu , Zhaohui Du PII: DOI: Reference:

S0045-7930(18)30701-1 https://doi.org/10.1016/j.compfluid.2018.09.021 CAF 4018

To appear in:

Computers and Fluids

Received date: Revised date: Accepted date:

30 January 2018 25 September 2018 28 September 2018

Please cite this article as: Chenxing Hu , Xiaojian Yang , Xiaocheng Zhu , Zhaohui Du , Stability and structural sensitivity analysis of the turbulent flow in the narrow vaneless diffuser with mean flow method, Computers and Fluids (2018), doi: https://doi.org/10.1016/j.compfluid.2018.09.021

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights The stability model is capable of predicting the frequency of instability.



The wavemaker region is localized at the centerline of the fully developed flow.



The existence of viscosity stabilizes the narrow vaneless diffuser flow.



The intersection region of the boundary layers may be effective for flow control.

AC

CE

PT

ED

M

AN US

CR IP T



1

ACCEPTED MANUSCRIPT

STABILITY AND STRUCTURAL SENSITIVITY ANALYSIS OF THE TURBULENT FLOW IN THE NARROW VANELESS DIFFUSER WITH MEAN FLOW METHOD Chenxing Hu, Xiaojian Yang, Xiaocheng Zhu, Zhaohui Du Shanghai Jiao Tong University, Shanghai 200240, China;

ABSTRACT

CR IP T

The performance and operating range of modern centrifugal compressors employed in the turbocharger are often limited by the vaneless diffusers. Despite of the simple geometry, the complex flow in the narrow vaneless diffusers makes it hard to identify the physicial mechanism of flow stability and requires more reliable model. In the present work, a three dimensional prediction model of the turbulent flow stability in the vaneless diffuser based on the mean flow method was proposed. The turbulent stability analysis was

AN US

performed by linearizing the incompressible Navier-Stokes equations around the numerically computed mean flows. The predictions of flow instability frequency and coherent structure were attempted. Then the adjoint stability analysis was carried out and the adjoint global mode was also obtained. Combining the direct and adjoint global modes, the structural sensitivity of the leading eigenvalue for the vaneless diffuser

M

was analyzed in order to localize the wavemaker region. Through validation against experimental data, the mean flow approach was proved to be capable of predicting the frequency of the flow instability in the

ED

narrow vaneless diffuser with mean flow appraoch. The receptivity analysis showed forcing or initial condition applied at the main flow upstream and the fully developed flow downstream the intersection

PT

region of the boundary layers on both walls of the vaneless diffuser may be most effective on the flow control. The structural sensitivity distribution indicated that the wavemaker region of the vaneless diffuser

CE

is localized at the centerline of the flow, especially near the radial position of 𝑟 =1.02~1.35. The viscosity

AC

and Reynolds stresses play a role of stabilizing the vaneless diffuser flow.

INTRODUCTION Centrifugal compressors play a significant role in some industrial applications and are usually

employed in turbochargers to boost internal combustion engines. The ever-increasing demand for more energy-saving compression system and components set new prerequisites for designers to meet the requirements without reducing the capacity of centrifugal compressors or excessively increasing the

2

ACCEPTED MANUSCRIPT

manufacturing costs. One major problem that requires special concern is the flow instabilities such as rotating stall and surge occurring at relatively low flow coefficient. Since the flow instability not only causes the deterioration in aerodynamic performance, but also are also the source of subsynchronous vibration and noise. The rotating stall in the vaneless diffuser, which often occurs before that in the impeller, may be the principal factor that limits the stable operation range of the centrifugal compressor in

CR IP T

the turbocharger.

Compared with vaned diffusers, vaneless diffusers usually have lower pressure recovery and wider stability margin. And vaneless diffusers are more suitable for off-design conditions and harder to come into choke than vaned diffusers [1-3]. In terms of instability, the physical mechanism of rotating stall in vaned

AN US

and vaneless diffusers are assumed to be different. For vaned diffusers, the vorticity shedding from the vane leading edge may convect back into the vaneless space and destabilize the vaned diffusers [4]. The nature of flow instability associated with radial vaneless diffusers have been extensively investigated by numerous researchers [5, 6]. For many efforts made by researchers to explain the mechanism and predict the occurrence of rotating stall within the vaneless diffuser, quantities of theoretical methodologies based on

M

different assumption and simplification have been proposed in the past decades. Generally speaking, one of the critical parameters may lie at the axial width of the vaneless diffuser. As shown in Fig.1, the influence

AC

CE

PT

ED

of the boundary layer and the main flow may differ significantly under narrow and wide vaneless diffuser.

(a) Narrow diffuser (b) Wide diffuser Figure 1. Sketch of flow structure in the vaneless diffuser Researchers such as Jansen [7, 8], Senoo [9, 10] and Frigne [11] assume that the three-dimensional wall boundary layer are supposed to be responsible for the occurrence of rotating stall in the narrow

3

ACCEPTED MANUSCRIPT

vaneless diffuser as shown in Fig 1 a). As is well known, the boundary layer on the vaneless diffuser walls is highly skewed owing to the pressure gradient and streamline curvature. The three-dimensional turbulent boundary layer may approach separation under certain inflow angle. Therefore, the researches based on boundary layer theory are mostly boundary layer calculation. And the calculations are associated with simplified boundary layer models and velocity profiles. In Senoo’s vaneless diffuser model, a boundary

CR IP T

layer calculation based on asymmetrical flow was conducted and two inflow angles were defined. An identified critical angle was defined for which radial reverse flow started in the vaneless diffuser and a critical angle is defined for which rotating stall occurred in the vaneless diffuser [9, 10]. Based on the numerical and experimental results, the critical angle corresponding to the rotating stall was 0.88 of the

AN US

identified critical angle. The critical angles for the vaneless diffusers with different axial width ratio were predicted with Senoo’s model. However, the prediction accuracy of the critical angle for narrow diffusers is relatively low based on the boundary layer theory. And this approach is highly depended on the choice of velocity profile and asymmetrical flow assumption. For the narrow diffuser with the fully developed flow, a turbulent stability analysis may be required.

M

Another theory about instability in vaneless diffuser is focused on the inviscid axisymmetric flow in the vaneless, which is adopted by Abdelhamind [12], Moore [13] and Shen [14]. As shown in Fig.1 b), the

ED

instability of inviscid main flow or the interaction between the inviscid main flow and boundary layers may be the principal reason for the flow instability in wide vaneless diffusers. In both works of Moore [13] and

PT

Shen [14], the diffuser flow instability were analyzed by fixing the coordinates in a frame rotating at the same speed as the stall cell, and a dense of resonant solutions were obtained where large pressure

CE

perturbations would occur regardless of the characteristics of the impeller. The inviscid model proposed by Hu [15] was realized by performing a linear stability analysis around inviscid flow. It was capable of

AC

predicting the occurrence of the stall in wide vaneless diffuser under different radius ratio, but failed to investigate the influence of the axial width of the diffuser. The large difference between point A and B as shown in Fig.2 requires more explanations. In summary, the limitations of inviscid models are that some phenomena and mechanism in physical sense cannot be clarified with inviscid analysis.

4

ACCEPTED MANUSCRIPT

16

A

12 10 8

Senoo's experiments Hu's Inviscid model

6 4

B

CR IP T

Critical Angle（ °）

14

2 0.00

0.04

0.08

0.12

0.16

Axial Width Figure 2. Stability margin of the vaneless diffuser of Rf=1.8 with Senoo’s experiments [9] and Hu’s inviscid model [15]

AN US

Given the two theories about diffuser stability, calculation of both the boundary layers and the main flow analytically in one time may lead to better understanding of the instability in vaneless diffuser. With the development of the linear stability theory, the stability method has been extended to turbulent flow such as cylinder wakes [16], airfoils [17] or confined jets [18]. For a typical turbulent stability analysis, the linearization of the Navier-Stokes equations around the mean flow is firstly conducted, and the dynamics of

M

the small amplitude perturbations around the mean flow is analyzed. The effect of the viscosity of the and the eddy viscosity

𝑡

ED

turbulent flow is generally reflected by adopting both the molecular viscosity

in the analysis. Based on the flow data obtained from experiment or numerical simulation by time averaging, which is called mean-flow approach, some researchers found that the nonlinear frequency and

PT

the structure of the coherent fluctuations can be predicted with linear stability analysis [19-22]. However,

CE

the work of Sipp and Levedev proved that these were not general phenomena and required certain conditions [23]. The work of Turton et al. indicated that the traveling wave case was successful but the standing wave case failed in predicting the nonlinear frequency with mean flow approach [24]. The

AC

physical relevance of the mean flow approach with the existence of a strong convective instability mechanism such as the K-H mechanism is revealed in the work of Beneddine. If a strong convective instability mechanism affects the mean flow at some frequency, the spatial structure of the modes should be very close to the spatial structure obtained [25]. For the vaneless diffuser concerned in the present work, the diffuser flow is associated with shear-layers and may be driven by Kelvin-Helmholtz instability [26]. Therefore, the mean flow approach for the vaneless diffuser flow is applied in the present work. Besides,

5

ACCEPTED MANUSCRIPT

the mean flow based on unsteady Reynolds-averaged Navier-Stokes simulation (URANS) can be applied for the turbulent stability prediction as indicated in Rukes’s work [27]. It should be noted that the choice of turbulence model has certain impact on the prediction accuracy. Since linear stability analysis usually provides the eigenvalues and the dominating dynamics of the system, another main problem we concern the most is how to locate the origin region of the instability

CR IP T

perturbations in the vaneless diffuser. In the past decade, the sensitivity study based on gradient method becomes a promising method in the hydrodynamics stability research and may provide more evidence of the instability mechanism. Adjoint operators, which is originally defined by Lagrange, has been thoroughly substantiated theoretically and broadly applied in solving many problems in mathematical physics. Under

AN US

the frame of local stability, adjoint method was first adopted as optimization techniques in fluid mechanics by Hill [28], and introduced for receptivity studies by Luchini [29]. By studying the Ginzburg–Landau model equation, Chomaz [30] showed that the wavemaker region could be identified as the overlapping region between the direct and adjoint global eigenvectors. While for the frame of global stability, the wavemaker region of the spiral vortex breakdown in swirling flows was the confirmed with structural

M

sensitivity analysis by Qadri and Juniper [31], and two different feedback mechanism related with the instability mode were revealed. The receptivity and structural sensitivity of global modes for cylinder wake

ED

was investigated with adjoint techniques by Gainnetti [32]. It was noticed that the wavemaker region was similar to the region where the placement of a control cylinder suppressed the vortex shedding in the

PT

experiments of Strykowski [33]. Sensitivity and receptivity study not only points out the wavemaker region of the flow field, but also may guide the flow control. In Fani’s work [34], both structural sensitivity and

CE

sensitivity to mean flow modifications can present quantitative indications on designing the control of the flow in symmetric channel with a sudden expansion.

AC

According to the above researches on stability and sensitivity, the quasi-laminar mixed method based

on eddy viscosity and molecular viscosity can be used to predict the instability of the shear flow and capture the dominating dynamics. Meanwhile, the structural sensitivity shows where an instability of the fluid flow is most sensitive to changes in the internal feedback mechanisms, and is capable of locating the wavemaker region of the vaneless diffuser. Applying the stability and sensitivity analysis to the vaneless diffuser flow could be beneficial on revealing the stability mechanism. The present study extends the

6

ACCEPTED MANUSCRIPT

previous inviscid core flow stability study of wide vaneless diffuser [15] and investigates the sensitivity of the turbulent flow in the vaneless diffuser. Main contents are organized as follows: The mean flow of two vaneless diffusers with different radius ratio is first sought by numerical simulation with two turbulence models, respectively. The direct global stability analysis is employed and the prediction results are validated against reported experimental data. Then an adjoint approach is carried out and the receptivity

CR IP T

based on the adjoint modes corresponding to the leading eigenvalues are investigated. Finally, the structural sensitivity analysis of instability mode in the vaneless diffuser is applied and the wavemaker region of the vaneless diffuser is localized.

FLOW CONFIGURATION AND THEORETICAL FORMULATION

AN US

Flow configuration and mean flow

The vaneless diffuser considered in this paper is a typical component in the centrifugal compressor used in turbochargers as shown in Fig. 3 a). In particular, an isolated vaneless diffuser with two parallel walls is considered. The impeller back sweep angle and the slip factor related to the impeller geometry is

M

omitted since they are beyond the range of the present study. A simplified configuration of the vaneless diffuser is reported in Fig.3 b). The highly distorted flow from the impeller flow into the vaneless diffuser

ED

with the radial velocity 𝑉𝑟 and circumferential velocity 𝑉𝜃 . An inflow circumferential angle 𝛼 is defined as 𝛼 = 𝑎𝑟𝑐𝑡𝑎𝑛(𝑉𝑟 /𝑉𝜃 ) at the diffuser inlet. When instability occurs in the vaneless diffuser, the

PT

corresponding inflow circumferential angle is defined as the critical angle 𝛼c . In what follows, the vaneless diffuser with a moderate radius ratio 𝑅𝑓 = 1.67 and the inlet Reynolds number 𝑅𝑒 =

AC

CE

105 is considered, where 𝑑1 is the diffuser inlet diameter, and 𝑉1 is the inlet velocity.

7

𝑉1 𝑑1⁄

= 3.6 ×

CR IP T

ACCEPTED MANUSCRIPT

(a) Compressor system (b) Vaneless diffuser Figure 3. Sketch of a vaneless diffuser in 3D model

AN US

Linearized N-S equations for turbulent flow

The fluid applied in the vaneless diffuser is assumed to be incompressible and described by the state vector 𝒒 = (𝑣𝑟 , 𝑣𝜃 , 𝑣 , 𝑝 ), where 𝑝 is the pressure and 𝒖 = (𝑣𝑟 , 𝑣𝜃 , 𝑣 ) are the radial, circumferential and axial velocity components, respectively. Then the non-dimensionalized 3-D, unsteady Navier-Stokes equations for the vanless diffuser flow are. ∂𝑡

+ 𝒖 ∙ ∇𝒖 = −∇𝑝 +

M

∂𝒖

1

∇2 𝒖

∇∙𝒖= 0

(1)

ED

(2)

Here, the variables are all non-dimensionalized according to Table 1. The velocity terms are non-

inlet radius 𝑟1 .

PT

dimensionalized by the impeller tip velocity U. The length is non-dimensionalized by the corresponding

CE

Table 1 Non-dimensional variables and parameters defined in terms of dimensional variables Definition

Nondimensional variable

Definition

Radius

𝑟 = 𝑟/𝑟1

Velocity

𝛎 = 𝛎/𝑈

Width

𝑧 = 𝑧/𝑟1

Pressure

𝑝 = 𝑝/(𝜌1 𝑈2 )

Time

𝑡 = 𝑡/(𝑟1 /𝑈)

AC

Nondimensional variable

For the fully saturated, turbulent shear flow we are dealing with, the flow is decomposed into a time-

̅ , a coherent perturbation or organized wave 𝒖 ̃ and a stochastic turbulent motion part 𝒖′ as averaged 𝒖 ̅+𝒖 ̃ + 𝒖′ given by Reynolds and Hussain [35]. In the present work, the interest is placed on the 𝒖=𝒖 ̃. coherent perturbations 𝒖

8

ACCEPTED MANUSCRIPT

Similar to the Reynolds-averaging procedure, the triple decomposition of velocity and pressure is substituted into the Navier-Stokes equations. Note that

𝜕̅ 𝜕𝜃

= 0 since the mean flow is circumferentially

symmetric. The governing equations for the mean flow field is as follows: 1

̅ ∙ ∇𝒖 ̅ = −∇𝑝 + 𝒖

′ 𝒖′ + ̅̅̅̅ ̅̅̅̅̅̅ ̅ − ∇ ∙ (𝒖 ̃⏟𝒖 ̃) ∇2 𝒖 ⏟ 𝒖

(3)

𝜏

𝜏

CR IP T

̅ = ∇ ∙ 𝒖′ = ∇ ∙ 𝒖 ̃=0 ∇∙𝒖

(4)

′ 𝒖′ + ̅̅̅̅ ̅̅̅̅̅̅ ̃𝒖 ̃ ) determines how the mean flow field is changed by the stochastic and coherent The term −∇ ∙ (𝒖 𝒖

Reynolds stresses. It can be seen that the mean-coherent and the mean turbulent interactions are already reflected in the mean flow field. The governing equations for the coherent perturbations are obtained as

∂𝒖 ̃ ∂𝑡

̃ ∙ ∇𝒖 ̅+𝒖 ̅ ∙ ∇𝒖 ̃ = −∇𝑝 + +𝒖

AN US

follows: 1

̃ + ∇ ∙ (𝒖 ̃𝒖 ̃ −𝒖 ̃𝒖 ̃ ) − ∇ ∙ (⏟〈𝒖′ 𝒖′ 〉 − ̅̅̅̅̅̅ ∇2 𝒖 𝒖′ 𝒖′ ) ⏟̅̅̅̅ 𝜏

(5)

𝜏

For the global linear stability analysis, emphasis on the small perturbation theory is made. It is indicated that any coherent disturbances of 2nd and higher orders are neglected, that is, the nonlinear term

M

̃𝒖 ̃−𝒖 ̃𝒖 ̃ is neglected. The remaining term 𝜏 represents the fluctuation of the stochastic turbulence 𝜏 𝑁 = ̅̅̅̅ 𝒖 due to the passage of the coherent wave. According to the work of Reynolds and Hussain [35], 𝜏 cannot be

ED

neglected and requires an appreciated model to obtain the closure equations. Hence, a Bossinesq’s turbulent hypothesis is adopted to invoke the closure. The time-mean and phase-average perturbations are shown in

CE

PT

equation (6):

2

〈𝒖′ 𝒖′ 〉 = 〈𝑘〉𝑰 − 2 3

𝑝 𝑡 (∇

+ ∇𝑇 )〈𝒖〉

2 ̅̅̅̅̅̅ ̅ 𝒖′ 𝒖′ = 𝑘̅𝑰 − 2 𝑡𝑡 (∇ + ∇𝑇 )𝒖

(6b)

3

with 𝑘 denotes the turbulence kinetic energy and

𝑡

(6a)

is the eddy viscosity. The superscripts 𝑡 and 𝑝

AC

denote a time-averaged and phase-averaged quantity. By time-averaging equation 6(a) and using the ̅ , equations (6) is transformed into 〈()〉 = () identity ̅̅̅̅̅ 2

2

̅ 𝜏̃ = 〈𝒖′ 𝒖′ 〉 − ̅̅̅̅̅̅ 𝒖′ 𝒖′ = 3 〈𝑘〉𝑰 − 3 𝑘̅𝑰 − 2 𝑡 (∇ + ∇𝑇 )〈𝒖〉 + 2 𝑡 (∇ + ∇𝑇 )𝒖

(7)

It is assumed that phase and time averaging changes the structure of the fine scale turbulence, but not ̅=𝒖 ̃ , equation (7) turns into the amplitude of the fluctuations. Therefore, 〈𝑘〉 = 𝑘̅. Considering 〈𝒖〉 − 𝒖 ̃ 𝜏̃ = −2 𝑡 (∇ + ∇𝑇 )𝒖

(8)

9

ACCEPTED MANUSCRIPT

After introducing equation (8) into the governing equations, the Reynolds number based on molecular viscosity needs to be redefined due to the introduction of eddy viscosity. An effective Reynolds number 𝜌𝑈𝐷

defined as 𝑅𝑒eff = 𝜇+𝜇 is employed here, and the final form of the linearized Euler’s equations for the quasi-laminar mixed approach are as follows: ∂𝑡

̃ ∙ ∇𝒖 ̅+𝒖 ̅ ∙ ∇𝒖 ̃ = −∇𝑝̃ + ∇ ∙ , +𝒖 ̃=0 ∇∙𝒖

1

(∇ + ∇𝑇 )𝒖 ̃-

(9)

CR IP T

∂𝒖 ̃

(10)

̃ = 𝟎 at the vaneless diffuser inlet and walls, and 𝑝̃ = 0 With the homogeneous boundary condition 𝒖 at the outlet. Generalized eigenvalue problem

AN US

Generally, the separation between spatial and temporal coordinates allows the use of Fourier modes in time when the variables are homogeneous. In the present analysis, the assumption of the circumferential uniform mean flow allows the coherent perturbation variables to be express in normal modes. Taking the quasi-laminar mixed approach for instance, the velocity and pressure perturbations are decomposed into the

M

form as follows:

(11)

𝑝̃ = 𝑝̂ (𝑟, 𝑧)𝑒 𝑖(−𝜔𝑡+𝑚𝜃)

(12)

ED

̃=𝒖 ̂ (𝑟, 𝑧)𝑒 𝑖(−𝜔𝑡+𝑚𝜃) 𝒖

̂ and 𝑝̂ are the amplitude function. 𝜔 is the complex frequency and 𝑚 is the Where 𝒖

PT

circumferential wave number. Since the complex frequency 𝜔 can be written as 𝜔 = 𝜔𝑟 + 𝑖𝜔𝑖 . If the imaginary part 𝜔𝑖 > 0, the disturbances will grow exponentially with time, and the mode will be unstable.

CE

If 𝜔𝑖 < 0, the disturbances will decay with time and the mode will be stable. And the rotating speed of the

AC

perturbation along the circumferential direction is related to the real part of 𝜔. 𝑓𝑟𝑠 =

𝜔

(13)

𝑚

After the norm form of perturbations in equations (11)-(12) are substituted into the equations (9)-(10),

the homogeneous ordinary differential equations can be expressed as 𝑴𝜙 = 𝜔𝑱𝜙

(14)

̂ , 𝑝̂ )𝑇 , 𝑴 and 𝑱 are listed in Annex A. Solving the equation (14) leads to a generalized where 𝜙 = (𝒖 eigenvalue problem, and the eigenvalue 𝜔 is used to determine the instability of the vaneless diffuser.

10

ACCEPTED MANUSCRIPT

Adjoint method With the calculated direct eigenmodes, the adjoint eigenmodes are essential to obtain the sensitivity of the eigenvalues to system internal or external changes. For considering the gradient-based sensitivity analysis with adjoint method, the adjoint quantities is defined as 𝒖† * 𝑝†

(15)

CR IP T

𝒒† = (

where 𝒖† and 𝑝† denote the adjoint velocity and pressure coherent perturbation.

Considering the axisymmetric nature of the flow configuration, the concept of standard 𝐿2 inner product in cylinder coordinate is introduced:

(𝒖, 𝒗) = ∫Ω 𝒖 ∙ 𝒗 𝑑𝑉 = 2𝜋 ∬ 𝒖 ∙ 𝒗𝑟𝑑𝑟𝑑𝑧

(16)

AN US

where centered dot refers to the usual scalar product of two vectors. The linear Euler’s operator, ℒ̃ has its corresponding adjoint operator ℒ † , which is defined in terms of the inner product as follows: (𝒒† , ℒ̃ 𝒒 ̃) = (ℒ † 𝒒† , 𝒒 ̃)

(17)

Then the adjoint governing equations can be derived from the following equations by integrating by

M

parts, which is called the continuous adjoint method.

𝒒𝑟𝑑𝑟𝑑𝑧 = ∬ ̃ 𝒒 ∙ ℒ † 𝒒† 𝑟𝑑𝑟𝑑𝑧 ∬ 𝒒† ∙ ℒ̃̃

(18)

ED

The boundary condition of the adjoint linearized N-S equations is obtained by ensuring the boundary terms during the integration all equal to zero. The detailed adjoint linearized N-S equations and boundary

PT

conditions are listed in Annex A.

With the similar treatment of the disturbance to the direct approach, the adjoint disturbance part of the

CE

variables can be written as

̂† 𝑒 𝑖(𝑚𝜃−𝜔 𝒒† = 𝒒

𝑡)

(19)

AC

̂† denotes to the adjoint modes. Then the eigenvalue problem for the calculation of adjoint stability where 𝒒 is derived as

̂† = 𝜔 † 𝑩† 𝒒 ̂† 𝑨† 𝒒

(20)

The same algorithm as the direct approach can be applied to calculate Eq. (20), and the adjoint eigenvalues and eigenmodes can be obtained. In order to value the adjoint eigenmodes directly, the adjoint ̂†, 𝒖 ̂ 〉 = 1. velocity modes are normalized with the condition 〈𝒖

11

ACCEPTED MANUSCRIPT

Receptivity and structural sensitivity The sensitivity of the direct global mode to structural changes in the linearized N-S equations is proportional to the overlap between the direct and adjoint global modes and defined as the structural sensitivity. The structural sensitivity mainly considers the effects on the system’s state of a variation in one of its own parameters, especially the disturbance field. The high structural sensitivity region, or the so

the largest drift of the eigenvalue.

CR IP T

called ―wave maker region‖ indicates the area where a modification in the structure of the system produces

Consider the linearized N-S equations with the eigenvalue 𝜔1 . And the source terms 𝛿𝑯 and 𝛿𝑅 are introduced to account for a possible physical forcing mechanism as shown in follows: ̂ ′′ + ℒ ′′ 𝒖 ̂ ′′ + ∇𝑝̂ ′′ = 𝛿𝑯 −𝑖𝜔1′′ 𝒖

AN US

̂ ′′ = 𝛿𝑅 ∇∙𝒖

(21) (22)

where double primes indicate quantities satisfying the perturbed equations. Assuming that the eigenvalue ̂ = *𝛿𝒖 ̂ , 𝛿𝑝̂ +. Then a simple expansion in 𝜔1 is drifted with 𝛿𝜔1 and the eigenmodes perturbation 𝛿𝒒 terms of the solution of the unperturbed problem.

(23)

𝑝̂ ′′ = 𝑝̂ + 𝛿𝑝̂

(24)

ED

M

̂ ′′ = 𝒖 ̂ + 𝛿𝒖 ̂ 𝒖

𝜔1′′ = 𝜔1 + 𝛿𝜔1

(25)

After inserting Eq. (23)-(25) into Eq. (21)-(22) and applying the Lagrange identity, the perturbation of

CE

PT

eigenvalue can be expressed with the adjoint modes over the flow domain 𝒟: 𝛿𝜔1 =

̂1† ∙ 𝛿𝑯 + 𝑝̂1† ∙ 𝛿𝑅]𝑟𝑑𝑆 ∫𝒟 [𝒖 ̂1† ∙ 𝒖 ̂1 𝑟𝑑𝑆 ∫𝒟 𝒖

(26)

Once the structural perturbations 𝛿𝑯 and 𝛿𝑅 are specified, the associated eigenvalue shift can be

AC

obtained. Here we consider the most significant structural sensitivity of the eigenvalue with respect to a localized feedback. The first kind of feedback is proportional to the velocity and affects the momentum ̂ and 𝛿𝑅 = 0, where 𝑪 is the equations. Then the structural perturbations are in the form of 𝛿𝑯 = 𝑪 ∙ 𝒖 2×2 matrix of the coupling coefficients. If the feedback is localized in space, 𝑪 can be expressed at the position (𝑟0 , 𝑧0 ): 1 𝑪(𝑥, 𝑦) = 𝛿(𝑟 − 𝑟0 , 𝑧 − 𝑧0 )𝑪0 𝑟

12

(27)

ACCEPTED MANUSCRIPT

where 𝑪0 is here a constant coefficient matrix and 𝛿(𝑟 − 𝑟0 , 𝑧 − 𝑧0 ) denotes the Kronecker delta function. Therefore the eigenvalue drift due to the localized feedback mechanism proportional to velocity can be derived with the Laplace transform: |∫𝒟 𝒖 ̂1† ∙ 𝑪 ∙ 𝒖 ̂ 𝑑𝑆| |∫𝒟 𝒖 ̂1† ∙ 𝒖 ̂1 𝑟𝑑𝑆|

≤ ‖𝑪0 ‖𝜆

(28)

CR IP T

|𝛿𝜔1 | =

where 𝜆 is the eigenvalue’s structural sensitivity to feedback at a given point: 𝜆=

‖𝒖 ̂1† ‖‖𝒖 ̂1 ‖ |∫𝒟 𝒖 ̂1† ∙ 𝒖 ̂1 𝑟𝑑𝑆|

(29)

AN US

Receptivity analysis mainly aims at the process how one environmental disturbance produces an instability wave in the flow. The external forcing such as one acoustic or one vorticity wave may be introduced into the flow as momentum forcing or mass injection, and the instability wave is induced. With adjoint eigenmodes in hand, the effects of the generic initial conditions and forcing terms on the timeasymptotic behavior of the solution of linearized Eulers’s equations can be derived.

M

Consider the Eq. (21) where the structural perturbation 𝛿𝑯 is the momentum forcing 𝑭 and initial

ED

condition 𝒖𝑖𝑛 , and 𝛿𝑅 is the mass injection 𝑚. Due to the linearity, the effect of the each term on the global mode can be studied separately once the direct and adjoint eigenmodes are determined. It can be seen that the effect of the forcing functions or mass injection is only associated with the adjoint variables,

PT

which stands for the Green’s function for the receptivity problem [32]. Therefore, the receptivity to the

CE

momentum forcing and initial conditions or boundary conditions is expressed by the adjoint velocity modulus |𝒖† |, and the receptivity to mass injection is expressed by the adjoint pressure modulus |𝑝† | in

AC

the present study.

MEAN FLOW AND NUMERICAL RESOLUTION Numerical implement of mean flow An accurate mean flow is a prerequisite for the reliability of the global stability analysis. Generally, the

mean flow for the vast majority of cases of industry interest can be determined by numerical means such as DNS or experimental method, which are high resource-consuming. Rather than calculating the flow field with DNS, the time average of 3D unsteady solution of the Reynolds averaged Navier-Stokes equations

13

ACCEPTED MANUSCRIPT

(URANS), which is also called mean flow approach, can be more cheaply obtained. The turbulent flow in vaneless diffuser can be analyzed with this method according to the previously researches [22, 27]. An isolated vaneless diffuser with structured gird is selected for the 3D URANS simulation of mean flow. The CFD software CFX is employed [36]. The time step for the unsteady simulation is ∆𝑡 = 1 × 10−4 and the calculation time is 0.8s, which is sufficient for the flow development. The time-averaged

CR IP T

velocity and eddy viscosity at the meridian plane are exported as the mean flow.

In centrifugal compressor, the impeller discharge a strong distortion in axial and tangential direction, resulting in asymmetry inlet conditions for the diffuser. Therefore, the inlet velocity distribution is assumed to be uniform in circumferential direction and the turbulence level of 𝑇𝑢 = 10% is applied at the vaneless

AN US

diffuser inlet. In addition, the inlet velocity distribution in the width direction is assumed to be uniform in order to eliminate the distraction of extra factors. The total pressure and total temperature are given as the inlet boundary condition, while the static pressure is given as the outlet boundary condition. As reported in Rukes’s work [27], the eddy viscosity calculated with the Reynolds Stress Model is more suitable in the stability study for strongly swirling jets. In the present paper, the 𝑘 − 𝜀 and Reynolds Stress Model are

M

employed respectively to confirm the turbulence model choice in the vaneless diffuser flow. For comparison, the geometric and flow parameters for the vaneless diffusers experimentally studied in

ED

Senoo’s work [37] are employed in the present paper as shown in Table 2. Considering the radius ratio and the width of the vaneless diffuser, they are both narrow vaneless diffusers and the boundary layer at the

PT

walls will merge in some radial position of the diffusers. According to the work of Megerle, the flow instability can be simulated with different turbulence models [38]. However, due to the influence of the

CE

boundary conditions and wall reflections in the isolated vaneless diffuser, the Koopman modes of the 3D

AC

URANS results computed through DMD will not be analyzed in the present paper. Table 2 Geometric and flow parameters of the mean flow Parameters Values bz 0.02 Rf 1.67,1.4 Inlet Re 3.6×105 Ma 0.1 Turbulence model 𝑘 − 𝜀, 𝐿𝑅𝑅 − 𝑅𝑆𝑀

Calculation of eddy viscosity

14

ACCEPTED MANUSCRIPT

The accurate prediction of the stability is related with the distribution and magnitude of the eddy viscosity. During the simulation of the mean flow, different turbulent model requires different strategy on calculating the eddy viscosity. For the eddy viscosity model such as 𝑘 − 𝜀 and Reynolds Stress model, the formation of eddy viscosity are based on the Bossinesq’s turbulent hypothesis as follows. 3

𝑡 𝑆𝑖𝑗

i, j, k=1,2,3

(30)

CR IP T

1 ̅̅̅̅̅̅̅ ′ ′ ′ ′ ̅̅̅̅̅̅ −𝑢 𝑖 𝑢𝑗 + 𝑢𝑘 𝑢𝑘 𝛿𝑖𝑗 =

However, the two-equation turbulence models such as 𝑘 − 𝜀 do not account for the transport of the turbulent shear stress, and the eddy viscosity is obtained based on the hypothesis of isotropic turbulent flow. The calculation of the eddy viscosity for the 𝑘 − 𝜀 turbulence model is completed as follow: 𝑡

= 𝐶𝜇 𝜌

𝑘

(31)

𝜀

AN US

where 𝑘 is the turbulence kinetic energy and 𝜀 is the turbulence dissipation rate, 𝐶𝜇 is constant.

While the Reynolds Stress Model, which solves both the normal and shear stresses, can be applied for the anisotropic turbulent flow. The calculation of the eddy viscosity requires a least-square solution of equation (30), and the eddy viscosity for the Reynolds Stress Model is =

(

M

𝑡

̅ ̅ ̅̅̅̅̅̅̅ (−𝑢 𝑢 + ̅̅̅̅̅̅̅̅ 𝑢 𝑢 𝛿 )∙( + ) ̅

+

̅

)∙(

̅

+

̅

)

𝑖, 𝑗, 𝑘, 𝑙, 𝑚 = 1,2,3

(32)

ED

Spatial discretization and stability calculations Once the mean flow has been obtained through the URANS method, the flow field data is then

PT

interpolated on a coarser mesh on which the three-dimensional generalized eigenvalue problems (14) are solved. The calculation domain of the vaneless diffuser are discretized in space by a spectral collocation

CE

method. 𝑁𝑟 + 1 and 𝑁 + 1 Chebyshev-Gauss-Labatto (CGL) points are chosen along the r and z direction. Then the generalized eigenvalue problem is calculated with QZ algorithm, which is realized by

AC

the function eig in the numerical software MATLAB. To eliminate the influence of grid number on the calculation, three different spatial resolutions are

employed to solve the eigenvalue problem with the same mean flow. In particular, the vaneless diffuser with the radius ratio of 1.67 at a large inflow angle is used for resolution verification. The leading eigenvalues (with the largest imaginary part) under three resolutions, which are 15 × 10(low-res), 35×25(mid-res) and 50×40(high-res) corresponding to radial and axial grid points, are shown in Table 3. It

15

ACCEPTED MANUSCRIPT

can be seen that the eigenvalues under middle and high resolution remains the same. Then the grid number of 35×25 is selected for the present calculation. Table 3 Leading eigenvalues calculations at three different spatial resolutions Case Eigenvalues Low-res 1.312-1.856i Mid-res 1.335-2.058i High-res 1.336-2.058i

CR IP T

STABILITY ANALYSIS Mean flow and eddy viscosity

According to the experimental results reported in Ref [37], the critical angles of the vaneless diffuser with the radius ratio of 1.67 and 1.4 are about 5.5° and 5.2°. In particular, when the inflow angle is 5.5°, the mean flow obtained by 𝑘 − 𝜀 model for the vaneless diffuser with the radius ratio of 1.67 and the axial

AN US

width of 0.02 is illustrated in Fig. 4. It can be seen that the angular momentum decrease gradually along the radial direction due to the strong viscosity effect near the wall. The circumferential velocity decrease from the inlet to the outlet due to the blockage effect of the skewed boundary layer.

In Fig. 4 (b), the distribution of the turbulence kinetic energy, which is calculated from the normal

M

Reynolds stress, is highly associated with development of the boundary layer on both shroud and hub. The intersection of both boundary layers occurs near the radial position of 𝑟 = 1.1. It can be seen that under the

ED

axial width of 0.02 which is quite narrow for the vaneless diffuser, the concept of the main flow as indicated by Moore is limited to the inlet of the vaneless diffuser. The turbulence kinetic energy is stronger

PT

near the wall rather than that at the centerline, and it decrease in the fully developed flow along the radial direction. While in Fig. 4 (c), the eddy viscosity comes to its maximum upstream the intersection near the

AC

CE

centerline of the diffuser. Then it decreases along the radial direction in the vaneless diffuser.

16

CR IP T

ACCEPTED MANUSCRIPT

AN US

(a) Angular momentum (b) Turbulence kinetic energy (c) Eddy viscosity Figure 4 Contours of the mean flow for the vaneless diffuser (Rf=1.67, bz=0.02) with 𝒌 − 𝜺 model

The mean flow calculated with Reynolds Stress model is shown in Fig. 5. Compared with that obtained with 𝑘 − 𝜀 turbulence model, the main difference of the mean flow calculated with Reynolds Stress model lies at the distribution of the eddy viscosity. Except for the eddy viscosity, the flow angle, turbulence kinetic energy and the angular momentum are much alike with those in Fig. 4. The reason for the difference

M

between the eddy viscosity of two model lies at that the shear stresses are considered in the calculation of the eddy viscosity with the Reynolds Stress model. The high eddy viscosity region is located at the

ED

centerline of almost the whole vaneless diffuser. This is consistent with the phenomenon that Reynolds

AC

CE

PT

stresses are anisotropic on the centerline in the classic channel flow.

(a) Angular momentum (b) Turbulence kinetic energy (c) Eddy viscosity Figure 5 Contours of the mean flow for the vaneless diffuser (Rf=1.67, bz=0.02) with 𝑹𝑺𝑴 model

17

ACCEPTED MANUSCRIPT

In order to further study the mean flow calculated with RSM model, the six Reynolds stresses calculated with the Reynolds transport equations are illustrated in Fig.6. It can be seen that the normal stresses are mostly distributed near the wall from inlet to the radial position 𝑟 =1.4. The normal and shear stresses are relatively small at the centerline of this region (𝑟 =1~1.4) compared with those near the walls. The magnitude of the stresses shows that the normal stresses in Fig.6 (a)-(c) are at least one order larger

CR IP T

than that of the shear stresses in Fig.6 (d)-(f). The shear stress may play little effect on the stability of the

M

AN US

vaneless diffuser.

ED

(a) 𝒗𝒓 ′ 𝒗𝒓 ′ (b) 𝒗𝜽 ′ 𝒗𝜽 ′ (c) 𝒗𝒛 ′ 𝒗𝒛 ′ (d) 𝒗𝒓 ′ 𝒗𝜽 ′ (e) 𝒗𝒓 ′ 𝒗𝒛 ′ (f) 𝒗𝜽 ′ 𝒗𝒛 ′ Figure 6. Contours of the Reynolds stresses for the vaneless diffuser with 𝑹𝑺𝑴 model Frequency and growth rate

PT

The stability characteristics of the mean flow are determined by the leading eigenvalue calculated by the generalized eigenvalue problem. The eigenspectrum of the vaneless diffuser with the radius ratio of

CE

1.67 is illustrated in Fig. 7. It can be seen that the flow field is still stable since the imaginary part of the eigenvalue is below zero. Compared with the continuous branches of Sipp and Lebedev [23], the imaginary part of the leading eigenvalue is a little away from zero. This is because that spectrum in Figure 7 is

AC

obtained with time-averaged simulation data under the critical angle from Senoo’s experiments, and the vaneless diffuser flow is still stable. The obtained spectrum can be seen as close to marginally stable considering the damping effect of the turbulence model. The stabilization effect of the viscosity may be over predicted with the present model.

18

CR IP T

ACCEPTED MANUSCRIPT

(a) 𝒌 − 𝜺 model (b) Reynolds Stress model Figure 7. Eigen-spectrum of the vaneless diffuser (Rf=1.67, bz=0.02) with inflow angle 5.5°

The predicted growth rate under different wave number 𝑚 of the vaneless diffuser with the radius

AN US

ratio of 1.67 and 1.4 are illustrated in Fig.8. It is assumed that the predicted azimuthal wavenumber of the leading eigenvalues refers to the stall cell number in the vaneless diffusers. It can be seen that the stall cell number are significantly influenced by the radius ratio. When the largest growth rate is obtained, the wave number of the vaneless diffuser with the radius ratio of 1.67 is 4. While for the vaneless diffuser with radius

M

ratio of 1.4, it’s 6. Considering the effect of the turbulence models on the prediction, the growth rate under different wave number with both turbulence model has the same trend. The reason why growth rate gives

ED

relevant results regarding the dominant azimuthal wavenumber is that the obtained growth rate in the present work is close to marginally stable, but may not be under the same strong conditions as the work of

PT

Sipp and Lebedev [23]. And the growth rate in the present work still carries the information for the

-0.6

Growth Rate

AC

CE

dominant azimuthal wavenumber.

-0.8 -1.0 -1.2 Rf=1.67 k-epsilon Rf=1.67 RSM Rf=1.4 k-epsilon Rf=1.4 RSM

-1.4 -1.6 -1.8

1

2

3 4 5 6 Wave Number

7

Figure 8. Growth rate under different wave number of both vaneless diffusers

19

ACCEPTED MANUSCRIPT

The predicted rotating speeds of the perturbation under different wave number are shown in Fig.9. Under small wave numbers, the rotating speed with negative values indicates that counter rotating stall may occur, which were also observed in experiments previously. The rotating speed of the diffuser with small radius ratio are generally a little larger than that of vaneless diffuser with large radius ratio. The reason for this may lie at that the angular momentum in the diffuser with large radius ratio suffers more loss due to the

CR IP T

increased friction loss. And the circumferential rotating speed of perturbation is highly related to the circumferential velocity, which contributes the most to the angular momentum. Besides, the influence of the turbulence model on the rotating speed of the leading eigenvalue are relatively small.

AN US

0.2 0.1 0.0

Rf=1.67 k-epsilon Rf=1.67 RSM Rf=1.4 k-epsilon Rf=1.4 RSM

-0.1 -0.2 -0.3

1

M

Rotating Speed

0.3

2

3 4 5 6 Wave Number

7

ED

Figure 9. Rotating speed of instability perturbation under different wave number Prediction results validation

PT

The predicted rotating speeds corresponding to the leading eigenvalues of the vaneless diffuser are validated against the experimental results in reference [31] in table 4. For comparison, the predicted results

CE

with the inviscid prediction model proposed previously are also included [14]. It can be seen that the predicted rotating speeds of the perturbations for both vaneless diffusers are close to those of the

AC

experimental results. The prediction accuracy of the mean flow obtained by 𝑘 − 𝜀 and Reynolds Stress models are close. Therefore, the following analysis is carried out based on the mean flow calculated with 𝑘 − 𝜀 turbulence model. Table 4. Leading eigenvalues validation against experimental results Radius ratio 1.67

Turbulence model 𝑘−𝜀 RSM Inviscid model

Re (𝜔) 0.17 0.162 0.435

Exp. results 0.167 0.167 0.167

Imag (𝜔) -0.81 -0.548 0.191

20

Reference value 0 0 0

ACCEPTED MANUSCRIPT

1.4

𝑘−𝜀 RSM Inviscid model

0.283 0.245 0.504

0.3 0.3 0.3

-1.049 -0.95 0.23

0 0 0

Direct modes analysis The eigenmodes of the leading eigenvalues indicates the dominant dynamics of the flow field and the structure of the perturbations. In order to reveal the stall mechanism in narrow vaneless diffusers, the

CR IP T

contours of the velocity perturbation corresponding to the leading eigenvalues of the vaneless diffuser are illustrated in Fig.10. In the term of the magnitude, the axial velocity perturbation is much smaller than the radial and circumferential velocity perturbations. As shown in Fig. 10 (a), a strong reverse perturbation from the outlet to the inlet exists in the fully developed flow field, which indicates that outlet reverse flow may exist before the occurrence of stall. The axial velocity perturbation has opposite propagation direction

AN US

upstream and downstream the intersection of the boundary layers. In Fig. 10 (c), the circumferential velocity perturbations may propagate in opposite directions and have the largest magnitude at the inlet and

PT

ED

M

outlet of the vaneless diffuser.

CE

(a) Radial velocity (b) Axial velocity (c) Circumferential velocity Figure 10. The direct eigenmodes corresponding the leading eigenvalue of the vaneless diffuser (Rf=1.67, bz=0.02)

AC

ADJOINT MODES AND SENSITIVITY ANALYSIS Adjoint modes and receptivity analysis Compared with direct global modes, the adjoint modes provide the gradient information of the vaneless

diffuser flow. The adjoint velocity perturbations corresponding to the leading eigenvalue are illustrated in Fig. 11. The strongest perturbation of the adjoint radial velocity is observed upstream and downstream the intersection region of boundary layer. Meanwhile, two adjoint circumferential velocity perturbations with

21

ACCEPTED MANUSCRIPT

opposite propagation directions are captured upstream and downstream the intersection. A comparison of Fig. 10 and 11 indicates a spatial separation of the direct and adjoint global modes related to the nonnormality of the vaneless diffuser flow. In addition, the intersection between the boundary layer on both

AN US

CR IP T

shroud and hub plays a significant role on both the direct and adjoint global modes.

(a) Radial velocity (b) Axial velocity (c) Circumferential velocity Figure 11. The adjoint eigenmodes corresponding the leading eigenvalue of the vaneless diffuser (Rf=1.67, bz=0.02)

M

Another important value of the adjoint mode is that the modulus of the adjoint velocity and pressure represent the Green function for the receptivity of the corresponding global mode. The scalar product of the

ED

adjoint mode with any forcing function or initial condition provides the amplitude of the direct mode [27]. Under this circumstance, the modulus of the adjoint velocity |𝑢† | and adjoint pressure |𝑝† | are shown in

PT

Fig. 12. It can be observed that the receptivity in both figures are symmetric and take effect near the inlet of the vanless diffuser. In many situations we are particularly interested in the effect of the initial condition or

CE

external forcing on the stability. In Fig. 12(a), the adjoint velocity indicates that the receptivity to initial condition or momentum forcing reaches its minimum at the intersection of the boundary layers. And it

AC

reaches its maximum at two regions, which are the centerline of the main flow upstream (𝑟 =1.07) and the fully developed downstream (𝑟 =1.11). Any momentum forcing or initial condition change applied on these two regions may be most effective on changing the eigenvalue of the system. While adjoint pressure in Fig. 12(b) shows that two regions in the flow near the intersection of the boundary layers have the largest receptivity to mass injection. Compared with the adjoint velocity, high adjoint pressure region is not limited at the centerline of the main flow upstream and fully developed flow

22

ACCEPTED MANUSCRIPT

downstream. The flow near the wall at the radial position of 𝑟 =1.02~1.08 and 1.15~1.27 is also receptive to mass injection. It indicates that when forcing such as flow suction or injection is applied on the high receptivity region of the vaneless diffuser, control of the flow stability may be obtained. And this provides an explanation for the effective passive control of the vaneless diffuser flow by introducing jet into the flow as reported in Tsurusaki’s work [39]. Meanwhile, the non-uniform distribution of the high receptivity

CR IP T

region along the width direction explains the different control effect when changing the injection position in

M

AN US

the width direction.

ED

(a) |𝒖† |

(b) |𝒑† |

Figure 12. The contours of the receptivity of the vaneless diffuser (Rf=1.67, bz=0.02) Structural sensitivity analysis

PT

With the direct and adjoint global mode, the wavemaker region of the instability perturbation can be obtained. The overlap region of the direct and adjoint global mode corresponding to the leading eigenvalues

CE

of both vaneless diffuser are illustrated in Fig. 13. From Eq. (34) and following Giannetti [32], it can be drawn that the core of the instability in both vanless diffusers is localized at the centerline of the turbulent

AC

flow, especially the region near the radial position of 𝑟 =1.02~1.35, which is near the intersection region of the boundary layers of the shroud and hub. Combining with the distribution of the Reynolds stress of the mean flow, it can be drawn that any local

feedback applied on the region with large normal Reynolds stress plays little effect on modifying the eigenvalue of the system according to the structural sensitivity distribution. In another word, the velocity

23

ACCEPTED MANUSCRIPT

distribution is reconstructed by the Reynolds stresses and viscosity near the wall, and the vaneless diffuser is stabilized. In addition, the structural sensitivity is relatively high at the centerline of fully developed region near the outlet for the vaneless diffuser with radius ratio of 1.67 compared with that in the vaneless diffusers with the radius ratio of 1.4. This indicates that the increased radius for the vaneless diffuser with the same

CR IP T

axial width brings higher possibility of the occurrence of self-excited perturbations. And this may be the

M

AN US

reason why the instability is easier to occur in the vaneless diffuser with larger radius ratio.

PT

CONCLUSIONS

ED

(a) Rf=1.67 (b) Rf=1.4 Figure 13. The contours of the structural sensitivity of both vaneless diffuser

With the turbulent mean flow, a stability and sensitivity analysis based on a quasi-laminar mixed

CE

approach is performed to compute the eigenmodes and sensitivities of the flow in the narrow vaneless diffusers. The prediction of rotating speeds of instability perturbations are validated against experimental

AC

results. By combining information from the direct and adjoint global modes, new physical insight for the occurrence of the instability in narrow diffusers is proposed. Some conclusions are drawn as followings: 1.

The proposed turbulent stability model based on quasi-laminar mixed approach is capable of predicting the rotating speeds of the instability perturbation in narrow vaneless diffusers.

2.

Considering the similar prediction accuracy of both turbulence models and the relatively small magnitude of shear stress, the isotropic eddy viscosity model such as 𝑘 − 𝜀 model is suitable for the calculation of vaneless diffuser mean flow.

24

ACCEPTED MANUSCRIPT

3.

Based on the adjoint modes, it can be drawn that two regions near the intersection of boundary layers, which are the centerline of the main flow upstream and the fully developed flow downstream, are most receptive to initial condition or external forces. The region near the wall at the radial position of 𝑟 =1.02~1.08 and 1.15~1.27 are also receptive to mass injection. The wavemaker region of the narrow diffuser with the axial width of 0.02 is localized at the centerline of the fully developed flow.

5.

CR IP T

4.

The existence of viscosity, especially the molecular viscosity and Reynolds stress, plays a stabilization effect on the stability of the narrow vaneless diffuser flow.

AC

CE

PT

ED

M

AN US

NOMENCLATURE r diffuser radius Rf ratio of outlet radius and inlet radius 𝑓𝑟𝑠 rotation speed of disturbances m circumferential wave number bz width ratio impeller blade geometry exit angle U impeller tip velocity 𝑣𝑟 radial velocity component 𝑣𝜃 circumferential velocity component 𝑝 static pressure ̅ 𝒖 mean flow velocity 𝑝̅ mean flow pressure ̃ 𝒖 velocity coherent perturbation 𝑝̃ pressure coherent perturbation † 𝒖 adjoint velocity coherent perturbation 𝑝† adjoint pressure coherent perturbation ′ 𝒖 velocity stochastic perturbation 𝑝′ pressure stochastic perturbation 𝛼 inflow circumferential angle for mean flow 𝑉𝑟 inlet radial velocity component for mean flow 𝑉𝜃 inlet circumferential velocity component for mean flow 𝛿𝑯 source term of external forcing applied on momentum equations 𝛿𝑅 source term of external forcing applied on conservation equation 𝑭 momentum forcing 𝒖𝑖𝑛 initial condition 𝑚 mass injection 𝑘 turbulence kinetic energy 𝜆 structural sensitivity † |𝒖 | receptivity to initial condition and momentum forcing

25

ACCEPTED MANUSCRIPT

|𝑝† |

receptivity to mass injection

REFERENCES [1] Deniz S, Greitzer M, Cumpsty A. Effects of Inlet Flow Field Conditions on the Performance of Centrifugal Compressor Diffusers: Part 2-Straight-Channel Diffuser. Journal of Turbomachinery,

CR IP T

2000, 122(1): V001T01A112. [2] Benini E, Tourlidakis. Design Optimization of Vaned Diffusers for Centrifugal Compressors Using Genetic Algorithms. AIAA Computational Fluid Dynamics Conference 2001, Anaheim, CA, 2001. [3] Japikse D. Centrifugal compressor design and performance. Concepts ETI, 1988.

[4] Everitt N, Spakovszky S. An Investigation of Stall Inception in Centrifugal Compressor Vaned

AN US

Diffuser. ASME 2011 Turbo Expo: Turbine Technical Conference and Exposition. 2012, 1737-1749. [5] Dazin A, Cavazzini G, Pavesi G, Dupont P, Coudert S. High-speed stereoscopic PIV study of rotating instabilities in a radial vaneless diffuser. Experiments in Fluids, 2011, 51(1):83-93. [6] Tsurusaki H, Munakata A. A study on the rotating stall in vaneless diffusers of centrifugal fans. 2nd report An analysis of pressure fluctuations. Nihon Kikai Gakkai Ronbunshu B Hen/transactions of the

M

Japan Society of Mechanical Engineers Part B, 1987, 53(487):885-893.

86(4): 750–758.

ED

[7] Jansen W. Rotating stall in a radial vaneless diffuser. ASME Journal of Basic Engineering, 1964,

607-617.

PT

[8] Jansen W. Steady fluid flow in a radial vaneless diffuser. Journal of Basic Engineering, 1964, 86(3),

CE

[9] Senoo Y, Kinoshita Y. Influence of inlet flow condition and geometries of a centrifugal vaneless diffuser on critical flow angle for reverse flow. Journal of Fluids Engineering, 1977, 99(1):98-103 [10] Senoo Y, Kinoshita Y. Limits of rotating stall and stall in vaneless diffuser of centrifugal compressors.

AC

ASME Paper 78-GT-19, 1978

[11] Frigne P, Braembussche D. A theoretical model for rotating stall in the vaneless diffuser of a centrifugal compressor, Journal of Engineer for Gas Turbines & Power, 1985, 107(2): 507-513.

[12] Abdelhamid N. Analysis of rotating stall in vaneless diffusers of centrifugal compressors. Canadian Aeronautics and Space Journal, 1980, 26: 118–128.

26

ACCEPTED MANUSCRIPT

[13] Moore K. Weak Rotating flow disturbances in a centrifugal compressor with a vaneless diffuser, Journal of Turbomachinery, 1989, 111: 442–449. [14] Chen H, Shen F, Zhu C, Zhao D. A 3D Compressible Flow Model for Weak Rotating Waves in Vaneless Diffusers—Part I: The Model and Mach Number Effects. Journal of Turbomachinery, 2011, 134(4): 041010.

CR IP T

[15] Hu C, Liu H, Zhu X, Chen H, Du Z. Dynamics of global instabilities in the vaneless diffuser: A numerical approach and its applications. Proceedings of the Institution of Mechanical Engineers Part G Journal of Aerospace Engineering, 2017:095441001771172.

[16] Pralits O, Brandt L, Giannetti F. Instability and sensitivity of the flow around a rotating circular

AN US

cylinder. Journal of Fluid Mechanics, 2010, 650: 513-536.

[17] Iorio C, González M, Ferrer E. Direct and adjoint global stability analysis of turbulent transonic flows over a NACA0012 profile. International Journal for Numerical Methods in Fluids, 2015, 76(3):147168.

[18] Juniper P. The effect of confinement on the stability of viscous planar jets and wakes. Journal of Fluid

M

Mechanics, 2010, 656(656):309-336.

[19] Meliga P, Pujals G, Éric Serre. Sensitivity of 2-D turbulent flow past a D-shaped cylinder using global

ED

stability. Physics of Fluids, 2012, 24(6):115. [20] Pier B. On the frequency selection of finite-amplitude vortex shedding in the cylinder wake. Journal

PT

of Fluid Mechanics, 2002, 458(458):407-417. [21] Barkley D. Linear analysis of the cylinder wake mean flow. Epl, 2006, 75(5): 194-205.

CE

[22] Oberleithner K, Paschereit O, Wygnanski I. On the impact of swirl on the growth of coherent structures. Journal of Fluid Mechanics, 2014, 741(2):156-199.

AC

[23] Sipp D, Lebedev A. Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows. Journal of Fluid Mechanics, 2007, 593(593): 333-358.

[24] Turton E, Tuckerman S, Barkley D. Prediction of frequencies in thermosolutal convection from mean flows. Physical Review E Statistical Nonlinear & Soft Matter Physics, 2015, 91(4):043009. [25] Beneddine S, Sipp D, Arnault A, Dandois J, Lesshafft L. Conditions for validity of mean flow stability analysis. Journal of Fluid Mechanics, 2016, 798:485-504.

27

ACCEPTED MANUSCRIPT

[26] Bousquet Y, Binder N, Dufour G, Carbonneau X, Trebinjac I, Roumeas M. Numerical Investigation of Kelvin–Helmholtz Instability in a Centrifugal Compressor Operating Near Stall. Journal of Turbomachinery, 2016,138(7): 071007-071007-9 [27] Rukes L, Paschereit O, Oberleithner K. An assessment of turbulence models for linear hydrodynamic

218.

CR IP T

stability analysis of strongly swirling jets. European Journal of Mechanics - B/Fluids, 2016, 59:205-

[28] Hill C. Adjoint systems and their role in the receptivity problem for boundary layers. Journal of Fluid Mechnics, 1995, 292: 183-204.

[29] Luchini P, Bottaro A. Gortler vortices: a backward-in-time approach to the receptivity problem.

AN US

Journal of Fluid Mechanics, 1998, 363:1-23.

[30] Chomaz M. Global instabilities in spatially developing flows: non-normality and nonlinearity. Fluid Mechanics, 2005, 2(37): 357-392.

[31] Qadri U, Mistry D, Juniper M. Structural sensitivity of spiral vortex breakdown. Journal of Fluid Mechanics, 2013, 720(4): 558-581.

M

[32] Giannetti F, Luchini P. Structural sensitivity of the first instability of the cylinder wake. Journal of Fluid Mechanics, 2007, 581(581): 167-197.

ED

[33] Strykowski J, Sreenivasan R. On the formation and suppression of vortex shedding at low Reynolds numbers. Journal of Fluid Mechanics, 1990, 218(218): 71-107.

PT

[34] Fani A, Camarri S, Salvetti V. Stability analysis and control of the flow in a symmetric channel with a sudden expansion. Physics of Fluids, 2012, 24(8): 084102.

CE

[35] Reynolds C, Hussain F. The mechanics of an organized wave in turbulent shear flow. Part 3. Theoretical models and comparisons with experiments. Journal of Fluid Mechanics, 1972, 54(2):263-

AC

288.

[36] ANSYS Inc. ANSYS CFX-Solver Theory Guide. 2006. [37] Kinoshita Y, Senoo Y. Rotating Stall Induced in Vaneless Diffusers of Very Low Specific Speed Centrifugal Blowers. Journal of Engineering for Gas Turbines & Power, 1985, 107(2):514-521. [38] Megerle B. Unsteady aerodynamics of low-pressure steam turbines operating under low volume flow conditions. Lausanne: EPFL, 2014.

28

ACCEPTED MANUSCRIPT

[39] Tsurusaki H, Kinoshita T. Flow Control of Rotating Stall in a Radial Vaneless Diffuser. Transactions

CR IP T

of the Japan Society of Mechanical Engineers Series B, 1999, 65(630):670-676.

ANNEX A

DETAILED FORMULATION OF THE EIGENVALUE PROBLEM

AN US

The generalized eigenvalue problem for the direct stability analysis is as follows: 𝑴𝜙 = 𝜔𝑱𝜙 𝑇 ̂ , 𝑝̂ ) , and the matrices 𝑴 and 𝑱 are listed as follows where 𝜙 = (𝒖 The matrix 𝑴 and 𝑱 are listed as follows 1 𝑑𝑣̅𝑟

𝑑𝑟

𝜕𝑟 ̅̅̅ 𝑣

+

𝑟

𝜕𝑟 𝑖𝑚𝑣 ̅̅̅

+

+

−𝜂

−

𝑟 1 2𝑖𝑚

𝑣̅𝑟

𝑟2

𝑅𝑒

2𝑣 ̅̅̅

𝜕 𝜕𝑟

𝑟

+

𝜕

𝑟

𝜕𝑧 𝑑𝑣̅𝑟

+

1

2𝑖𝑚 𝑟2

𝑅𝑒

𝑖𝑚𝑣 ̅̅̅

+

𝑟

𝑑𝑣̅𝑧

(

𝑖𝑚

𝑣̅𝑟 𝑟

−𝜂

𝑣̅𝑟

0

𝑑𝑟

M

𝑴=

𝜕

+ 𝑣̅𝑟

𝑑𝑟 𝑑𝑣 ̅̅̅

𝜕

+

𝑟

𝜕

𝜕𝑟

+

𝑖𝑚𝑣 ̅̅̅ 𝑟

+

(33)

0 𝜕

𝑑𝑧 𝑑𝑣 ̅̅̅

𝑑𝑣̅𝑧 𝑑𝑧

𝜕𝑟 𝑖𝑚

𝑑𝑧

+ 𝑣̅𝑧

𝜕 𝜕𝑧

−𝜂−

1

1

𝑅𝑒

𝑟2

𝑟 𝜕

𝜕𝑧 )

(34)

PT

ED

0 0 0 0 𝑖 0 0 0) 𝑱=( 0 𝑖 0 0 0 0 𝑖 0

1

.

𝜕

𝜕𝑟

−

𝑚

𝜕

𝑟

𝜕𝜃

CE

𝜂=

+

𝜕

+

𝜕

1 𝜕 𝑟 𝜕𝑟

−

1 𝑟

(35)

𝜕

/+

𝜕 𝜕𝑟

𝜕𝑟

+

𝜕

𝜕 𝜕

(36)

𝜕

AC

The adjoint linearized N-S equations are listed as follows:

𝜕𝑣

𝜕𝑡

+ .−

𝜕𝑣 𝜕𝑡

̅̅̅ 𝑣 𝑟

−

𝜕𝑣 ̅̅̅ 𝜕

𝑣 𝑟

−

1

1

/ 𝑣𝑟 † + .

𝑟 1

2 𝜕𝑣 𝑟

+ .−

2𝑣 ̅̅̅̅ 𝑟

/ 𝑣𝑟 † + .−

/

𝑟𝜕𝜃

𝜕𝑣 ̅̅̅ 𝜕

− 𝑣̅

+

𝜕𝑣

̅̅̅̅ 𝑣

+

𝑟

𝜕𝑟 𝜕𝑣 𝜕

𝜕𝑟

𝜕

+

+

𝜕𝑣 ̅̅̅̅

𝜕𝑣

− 𝑣̅ 𝜕𝑣 ̅̅̅

𝜕𝑟

1

𝜕𝑝

+

𝑟 𝜕𝜃

/ 𝑣𝜃 † + 1

− 1 𝑟

1 𝜕𝑣

.

𝜕𝑣

𝜕𝑣 ̅̅̅ 𝜕𝑟

𝜕 𝑣 𝜕𝑟

=0

𝜕

𝑣

†

1

.

𝜕 𝑣 𝜕𝑟

1

+ .−𝑣̅𝑟 + 1 𝜕 𝑣

+𝑟

𝜕𝜃

/ 𝑣𝜃 † + .−𝑣̅𝑟 +

− 𝑟𝜕𝜃 −

(37)

1

𝜕𝜃

29

𝑟

𝜕 𝑣 𝜕

1 𝜕𝑣 𝑟

1 𝜕 𝑣

+𝑟

+ /

+

1 𝜕𝑣

𝜕𝑟 𝜕 𝑣 𝜕

/

−

𝜕𝑟

𝜕𝑝 𝜕𝑟

+ .−𝑣 ̅̅̅𝜃 −

/=0 −

1

(38) 2 𝜕𝑣 𝑟 𝑟𝜕𝜃

/=0

−𝑣 ̅̅̅ 𝜃

𝜕𝑣 𝑟𝜕𝜃

−

(39)

ACCEPTED MANUSCRIPT

𝜕𝑣 𝜕𝑡

+

𝜕𝑣 ̅̅̅ 𝜕

𝑣𝑟 † +

𝜕𝑣 ̅̅̅̅ 𝜕

𝑣𝜃 † +

𝜕𝑣 ̅̅̅ 𝜕

𝑣

†

+ .−𝑣̅𝑟 +

1

1 𝜕𝑣 𝑟

/

𝜕𝑟

− 𝑣̅

𝜕𝑣 𝑟𝜕𝜃

−

𝜕𝑝 𝜕

−

1

.

𝜕 𝑣

1 𝜕 𝑣

+𝑟

𝜕𝑟

𝜕 𝑣

/=0

𝜕

+

𝜕𝜃

(40)

with the boundary condition Inlet: 𝑣𝑟 † = 0

†

= 0，𝑝† +𝑣𝑟 † (𝑣̅𝑟 −

= 𝑣𝑟 † = 𝑣𝜃 † = 0

AC

CE

PT

ED

M

AN US

Walls: 𝑣

†

1

1

) = 0 (41b)

CR IP T

Outlet: 𝑣𝜃 † = 𝑣

(41a)

30

𝑓

(41c)