Prediction of Stall Phenomenon in Blade-to-Blade Passage of Axial Flow Compressor

The flow field between two blades was analyzed numerically by solving the steady, two dimensional and incompressible (time-mass averaged Navier- Stokes equations). The ( κ - ε ) t urbulent model was used to simulate the condition of the momentum equation and to obtain the eddy viscosity. The SIMPLE algorithm is used to satisfy velocity-pressure coupling method and to satisfy the conservation of mass. A computer code was constructed in this work using FORTRAN 90 language to simulate the flow. The numerical results were compared with experimental results of other researcher for flow through the cascade of axial compressor, and were found to be in a good agreement. Three different airfoils common, used for axial flow compressors blades, were investigated in the present study. The study shows that stall happens. At incidence angles (-8 0 to 7 0 ) for NACA 65(12)10 and (-9 0 to 3 0 ) for NACA 65(18)10 and (-7 0 to 9 0 ) for NACA 65(8)10. The result also show that stall happens when the total pressure loss is greater than 0.06 for NACA 65(8)10 when stagger angle is greater than 45 0 stall happens for total pressure loss less 0.06 as the outlet flow angle decreases. However, the quasi-three dimensional, steady, incompressible, turbulent, adiabatic and single-phase fluid flow inside the blade-to-blade passage of an axial flow compressor stator was also studied and the study shows that Stall was seen to take place at blade tips during starting.


Introduction
Stability in a compressor is the ability of a compressor to recover from disturbances that alter the compressor operation about an operational equilibrium point. Disturbances may be considered as transient or deliberate changes to the operating point. A rapid increase in pressure across the blade causes a marked thickening of the boundary layers and produces an effective contraction in the flow, thus a contraction coefficient is introduced in the model line. In the case of transient disturbances, the system is stable if it returns to its original operating point. If the disturbances drive the compressor away from the original point, the system will be unstable. The steady state match between a compressor and its drive turbine or jet nozzle, which is perturbed by a transient change of mass-flow, is a good example of this case. When there are deliberate changes to the operating point, the performance is considered stable if a new operational equilibrium point can be achieved, e.g., shifting the operating point by changing the compressor shaft speed. If steady state operation at a new operating point is not possible, the system is unstable. Stability in compressors may be studied from two different perspectives. The first is called operational stability, which deals with the matching of compressor performance with a downstream flow device such as a turbine or throttle. The second is aerodynamic stability, which deals with deteriorations in the operation due to flow separation, stall or surge [1].

Stall
During the normal operation of a compressor, the airflow through the compressor is essentially steady and axisymmetric in a rotating coordinate system. If a flow, instability is somehow introduced into the system (say, due to a change in the rotor speed, flow separation at the inlet, or other type of flow distortion), instabilities may develop and the compressor performance may deteriorate. The instability manifests itself as either a rotating stall or surge. It often takes only a few seconds for rotating stall to build up and the compressor can operate under rotating stall for several minutes before damage develops. Rotating stall can occur in both compressible and incompressible flow. In a coordinate system attached to the blades, rotating stall moves in a direction opposite to the blade motion at a fraction of the rotor speed. However, in the inertial coordinate system, the stall region propagates in the same direction as the wheel motion. The number of stall cells depends on the compressor at hand; one to nine stalled cells has been reported. Two types of stall associated with the number of stalled cells exist, progressive and abrupt. In progressive stall, a phenomenon involving multiple stalled cells, the pressure ratio after stall reduces gradually. Abrupt stall results in a sudden drop in total-to-total pressure rise, and appears to always involve a single stalled cell. One of the characteristics of pure rotating stall is that the average flow is steady with respect to time, but the flow has a circumferentially non-uniform mass deficit. During rotating stall, the cyclical variation of the pressures on the blades can cause them to fatigue and eventually break. Several types of rotating stall exist [2]: • Part-Span: A restricted region of the blade passage, usually the tip is stalled. Stall near the root has also been reported.
• Full-Span: The entire height of the annulus is stalled.
• Small/Large scale: In this case, a small/large part of annular flow path is blocked.

Mathematical formulation:-Governing equations
The two-dimensional instantaneous governing equations of mass momentum for steady, turbulent, incompressible flow in a coordinates system can be written in tensor conservation form expressed in Cartesian coordinates as follows [3].
Where, t ij is the viscous shear stress tensor that is expressed as [4]: The Reynolds stress tensor τ ij can be determined according to the boussineq assumption as: Where : µ t : is the turbulence eddy viscosity and estimated by the ε − k two equation turbulence model : The differential form of the turbulence kinetic energy ( k ) and dissipation rate of turbulence kinetic energy ( ε ) have given as [5]: is the turbulence production term. The last coefficients appearing in equation (7) are as those adopted in ref. [4] in the standard k-ε two-equation turbulence model. These coefficients are:

Boundary conditions
At inlet, Cartesian velocity component (u, v) is prescribed with inlet flow angle. For turbulence quantities, such as ( k ) and ( ε ) are normally not known, but they must be estimated. Usually [6] At exit, the velocity distribution is decided by what is happing within the domain [7]. The gradients normal to the outlet surface of all quantities is assumed zero.
The wall is the most common boundary encountered in confined fluid flow problems. In this section, a solid wall parallel to the ξ-direction is considered. The impermeable no-slip condition ( v=u=0) is the appropriate condition for the velocity components at solid walls. The turbulence scalar transport equations for ( ε and k ) are only valid for fully turbulent regions. An additional model must be introduced to treat the laminar sublayer region. The (Wall Function Method) [4] is used in the present study to eliminate the large number of grid points needed to resolve the laminar sub-layer at near-wall regions. The following function are used to bridge the near -wall region: Here, y is the conventional coordinate normal to the wall and w τ is the wall shear stress. The subscript p refers to the grid node next to the wall; k and E are the constants from the law of the wall, with values of 0.4 and 8.8, respectively [4].

Transformation of the governing equations
The governing equations can be expressed in the general form [7] ( ) ( ) Eq. (12) changes according to general transformation ξ= ξ(x,y), η=η(x,y). partial derivative of any function f are transformed according to Where J is the jacobain of the transformation given by Since the strong conservative form of the equation is desirable for numerical computations, the integral form of the conservational equations over a finite volume element is preferable. Upon introducing The following integral conservative relation is obtained from eq. (12) for an arbitrary scalar dependent variable Φ:

Numerical formulation:-Discretization of the governing equations
The general form of the governing equation (17) is integrated over each discrete CV in the computational plane (ξ, η). Let a new working variable be (Iø) j such that superscript j can be any of the computational directions (j=ξ or η). This is called the total flux in the j th direction. It is expressed mathematically as follows: (17) gives:- For convenience, the above notation can further be simplified to:- The requirement now is to express the (I 's) at node (P) in equation (21) in terms of the properties of the considered (P) node and its neighbor nodes (E, W, N, and S).
As a typical analysis, the I e can be considered as follows: - Assume (F e ) is the mass convective flux at face (e), which can be defined as follows: Expanding the second term in equation (22) (by expressing the partial derivative w.r.t. ξ using central differencing scheme) and substituting for Γ 1 utilizing its definition given by the relation in equation (18) gives: -Blade Passage of Axial Flow Compressor At this point, it would be useful to define a quantity D e ; it is the diffusion term coefficient at face (e). Hence, this coefficient can be expressed as: - Substituting equation (25) in equation (24) and the result in equation (22) gives: - (26) contains the quantity ø e , which needs to be specified in terms of neighbor nodal values. This is usually achieved by using an appropriate interpolation scheme. In the present work hybrid scheme is used. for simplifying the analysis discussed the upwind scheme and after some arrangement is converted to hybrid scheme based on Peclet number.
The next step is that, the continuity equation is integrated over the same CV. (32) is multiplied by ø P to get: ( ) Hence, substituting the result of (31) in to (34) and rearranging gives: Where the Discretizing coefficient of node P can be define as follows: (36) is a general discretized form of the conservative general transport equation and can be applied for both a staggered and nonstaggered CV by setting the appropriate boundary at the interfaces. The argument of ø can hold for any dependent variables.

Grid arrangement for the dependent variables
The use of staggered grid in the body fitted coordinate system flows, requires the solution of momentum equations for covariant velocity components on the faces of the control volume. Calculating of coefficients and geometrical factors (a 1 , b 1 , a 2 ) must be performed for many sets of grid, hence programming becomes very difficult and more effort is required to calculate the coefficients and geometrical factors. These problems can be avoided through the solution of the momentum equations by using Cartesian velocity components on a non-staggered (collocated) grid where all variables are stored in the center of the control volume as in [8].

Pressure correction equation (P. C. E.):
After solving the momentum equation, the velocity field obtained does not guarantee the conservation of mass unless the pressure field is corrected. Therefore, the velocity components (u,v) ,(U,V) and the pressure must be corrected according to the continuity equation as follows:- Where superscript ( * ) refers to the last iteration values and the superscript (') refers to the corrected variables. the continuity equation is rewritten after introducing the correction of contravariant velocity components (U, and V in equation (37)) to get:- Equation (39) represents the discretized continuity equation called full pressurecorrection equation.
The cross-derivatives of pressure correction in equation (39), such as ( ne P′ , se P ′ , nw P′ , sw P′ ) are equal to zero for orthogonal grid. for non-orthogonal gird, if the cross-derivatives of the pressure correction are retained, the derived pressure correction equation will be a nine-diagonal matrix for twodimensional flows. Solving this matrix is complex and expensive, therefore the cross-derivatives terms are usually omitted for simplicity. This is called simplified pressure correction equation [9] and [10]. in the present work the crossderivative is calculated as [3] The u-Cartesian velocity at east face is usually obtained by liner interpolation, i.e.
In a collocated grid arrangement this may lead to pressure oscillations. To avoid this, the face velocities are obtained from (Rhie and Chow interpolation) [11] by subtracting and adding the pressure gradient to equation (43) as follows: The first term is calculated as a weighted average. The second term is represents a fourth-order derivative term to dampen oscillations. The v Cartesian velocity components at face (e) are simply averaged weighted linear interpolation without the need to any further treatment. This is because the pressure gradient across face (e) does not affect these components, hence:- In the same way can obtain the Cartesian velocity component for the control volume faces (w, n, s)

Pressure-velocity coupling method
In the present work, the (SIMPLE) algorithm (semi-implicit method for pressure-linked equations) is used to couple the pressure and velocity as in [12]. This method is done by solving the momentum equations using the guessed pressure field to obtain the velocity field. The velocity filed obtained satisfies the momentum equations. Then the velocity and pressure are corrected because the velocity field violates the conservation of mass.

Convergence criteria
At the end of each solver iteration, the residual sum for each of the conserved variables is computed and stored, thus recording the convergence history. The residuals decay to some small value and then stop changing. In this work an iteratively converged solution is assumed to be reached when the largest residual of all variables is less than 10 -4 , where residual can be define as follow:

Grid generation
The simplest grid generation technique is the algebraic method. The derivatives of the boundary in the physical plane provide even more flexibility in the mapping. For instance, orthogonally at the boundary can be forced in the physical plane. The interior grid points of the algebraic gird in the physical domain can be calculated by using the following algebraic equations: - Where (x (i, j), y (i, j)) are any interior points in the physical domain, m and n are the node numbers in the x and y directions respectively, and l is the axial length of the physical domain. The grid that generated algebraically by using equations (47) for the NACA 65(12)10 linear cascade is shown in Fig.  (1).

Quasi-three dimensional flow:-
The through-flow field will be applied by using quasi-three dimensional flow. The mathematical approach for describing quasi-three dimensional flow will be discussed using the radial equilibrium theory.

Radial equilibrium theory
The earliest approaches to radial flow analysis were based on the radial equilibrium theory.
The basic assumption of the radial equilibrium type of design is that the radial velocity (vr ) is zero at entry and exit from a blade row. The equation of motion using cylindrical coordinates for incompressible, inviscid flow is [13]: If there is a radial equilibrium, equation (48) may be written as: The total pressure P t is given by:  (52) will hold for the flow between blade rows of the compressor, and may be used to determine the axial velocity (u) variation once the tangential velocity (v) distribution is chosen. If (v = f(r)) then the distribution of (u) with radius is obtained.
The free vortex design satisfies the requirement of radial equilibrium [14], from this case obtains: r. v = constant and no variation of (u) across blade height.

Results and discussions
The computational model was built based on the linear cascade of NACA 65-(12) 10 compressor blade. In addition, to study the behavior of the compressor before and after stall has happen; two cases were taken. The treatments of boundary conditions of these cases were described in article (2.2) and the input data needed were recorded in table (1).
Case (1), the velocity vectors and the contours of velocity and static pressure at mid-span (at 50% of the blade's height) at best incidence angle are presented in Figures. (2), (3) and (4) respectively. Velocity vectors show that the fluid is changing its direction from inlet to exit following the signature of blade profile.
From figure (3) due to the shape of the airfoil, the flow is seen to accelerate along the suction and pressure sides near the leading edge, and the maximum velocity occurs at a position just off the suction side at about 21% of blade chord. These phenomenon prove that the compressor works as a nozzle-diffuser.
The static pressure contour figure (4) shows that, on the suction side, the pressure falls very rapidly from the stagnation point towards the throat, reaching a minimum value at 26% of chord line. After that acceleration occurs until reaching the trailing edge. On the pressure side, the same behavior can be seen but the position of the minimum is at 21% of chord line.
For case (2) figures (5), (6) and (7) represent the velocity vector and contours of the velocity and static pressure at mid-span when the stall happens. From figures (5)and (6) it is noticed that the direction of the velocity changes according to the airfoil of blade, but on the suction side the velocity is reduced to reach a minimum value (about 4 m/s) at 37% of chord line, then the velocity reverses its direction to form a vortex. Figure (7) Shows that the pressure decreases towards the exit direction. This is on the suction surface. On the pressure surface the pressure decreases until reaching minimum value at 42% of chord line. After that acceleration occurs until the trailing edge.
The incidence at which stalling occurs is difficult to define precisely, and the stalling point is usually arbitrary specified as the incidence at which the total pressure loss is less than or equal to 0.06. [14]. The incidence angle is plotted with total pressure loss as in figure (8). It is noticed that the total pressure loss after 0.06 increases rabidly and the minimum total pressure loss is at incidence of -1°.
The stall limit of NACA 65(12)10 for multi stagger angle, is illustrated in Figure (9), it is noticed from this figure that the range of inlet flow angle become little if stagger increases, and it is noticed that the outlet flow angle do not change at the same stagger angle. This is the same for (Carter rule) [14].
Figure (10) shows the relation between total pressure loss and diffusion factor. It is noticed that if the diffusion factor is larger than 0.6 the total pressure loss increases rapidly. So the stall limit is at diffusion factor 0.6 or less. But for larger camber angle NACA 65(18)10 the diffusion may be at over 0.6 and total pressure loss is less than 0.06 without causing stall to happen as shown in figure (11).
To explore the effect of the camber angle on the stall limits, the stall limit is plotted for NACA 65(18)10 and NACA 65(8)10 in figures (12) and (13) respectively. It is noticed that stall limit acts at different incidence angles.
Long compressor blades are twisted from hub to tip. There are many methods to calculate this twist. In the present work the free vortex theory was used. Boundary conditions for this problem are the same as those of case one. The hub to tip ratio taken is (0.5). The dimensions of the blade were, the chord length (15cm) and span length (30cm).
The inlet and outlet flow angles were plotted with span from hub to tip in figure (14). It is noticed that both inlet and outlet flow angles decrease from hub to tip but outlet flow angle has a large slop. Figure (15) presents the variation of the total pressure loss coefficient with span from hub to tip. It is noticed that the maximum total pressure loss act at the tip of the span. This means that the stall began at the tip of the blade.

The
Comparison between the experimental and theoretical results: The theoretical results obtained in this study by using numerical technique were compared with the experimental results of [sabah 2008]. Figure ( This figure shows that at the section 0.125 from the space, the present result lower than that of the reference. However, the same behavior is noticed. At the middle of the blade to blade spacing, the results are vary close. While, at 0.875 from space, the results are also close. This figure shows that the results are different at the suction surface and get closer as we approach the pressure surface. However, the behavior in all cases is nearly the same.
This comparison agreement is considered to be good bearing in mind, the approximation of theoretical results, as compared with actual practical results.