Abstract
The object of this paper is to provide a reliable tool to simulate the stationary and transient operation performances in multistage axial compression systems, especially poststall behavior. An adapted version of the 1D Euler equations with additional source terms is solved by a time marching and control volume method. The equations are discretized at midspan both inside the blade rows and the nonbladed regions, along the real flow path geometry. The source terms express the blade-flow interactions and are estimated by calculating the velocity triangles for each blade row. Loss coefficient and deviation models are supplied by empirical correlations and are compared to experimental data in all flow regions. Transient simulations are carried and compared to the experimental results for several values of parameter B. The flow mechanism inside the compressor during a surge cycle is also shown.
1. Introduction
The compression system is an important component of a gas turbine engine. During the operation of axial multistage compressors, instabilities such as surge and rotating stall were observed in many experimental tests. However, such tests are expensive and not practical for preliminary design purpose. Numerical simulations of poststall performances of compression system are needed and useful in the design phase. In this paper, a reliable numerical tool is presented to simulate the instability behavior and to analyze the flow mechanism inside the compressor during poststall process.
This tool is only based on the compressor geometry and empirical correlations to simulate the poststall performance. The mesh cells are not only distributed between blade rows, but also inside of each blade row, allowing for a detailed representation of the flow field. An adapted version of Euler equations, mass, meridional momentum, circumferential momentum, and energy balances being applied for each control volume at mean radius. These equations are solved by a time-marching, finite-volume method [1]. The source terms expressing the blade-flow interactions are estimated by calculating the velocity triangles for each blade row with the help of empirical correlations.
The test compressor rig of Day is chosen to validate this model because of the detailed experiment results [2]. The steady characteristic curve of the compressor is predicted and is compared to the experiments. Dynamic simulations of the compression system are also carried out to predict the poststall performance and to validate the present tool. The flow mechanism inside of the compression system during one surge cycle is studied and the frequency of the surge is analyzed.
2. Background
Rotating stall and surge are two kinds of poststall phenomena in compression systems [3]. Rotating stall is a severe nonaxisymmetric flow which includes one or several stall cells, propagating in the circumferential direction at a fraction of the rotor speed. Rotating stall can lead to serious mass flow and pressure ratio reductions. However, the annulus average mass flow and pressure ratio remain almost constant once the stall cells are fully developed. Surge is an unsteady 1D phenomenon with mass flow oscillations in the axial direction, even the reversed flow inside the compressor at severe conditions. Both rotating stall and surge can lead destructive to consequences to the compressor; however the appearance of the nonrecovery rotating stall in the aero-engine system, make it important to study the poststall behavior during the transient operations.
If only the global behavior is concerned, a lumped parameter technique can be implemented, which considers each component in the compression system as one control volume on which mass, momentum and energy balances are applied. The interactions of the compressor and the transient flow are included through source terms which are derived from steady characteristics. This method can provide fast and reliable results, but it lacks the ability to simulate the flow details inside the compressor. Greitzer developed such a lumped-volume method to study surge and rotating stall phenomena, based on the compressor steady characteristics [4, 5]. Davis applied a lumped volume method for each stage of the compressor and carried out parametric studies [6]. Morini developed a nonlinear modular dynamic model to study the deep surge of compression systems through an imaginary steady-state performance map [7].
If the mechanisms of stall inception and rotating stall development are the main focus, 2D and 3D models are necessary. Takata and Nagano developed a 2D model to simulate rotating stall in axial compressors, with the assumption that blade rows are replaced by semiactuator disks [8]. Moore and Greitzer proposed a 2D, unsteady, lumped parameter model to study the dynamic response of compression systems and the coupling between rotating stall and surge [9]. Lindau and Owen developed a quasi-3D actuator-disk model to simulate rotating stall and surge, with source terms representing blade losses and deviations [10]. Escuret used 3D Euler equations with source terms estimated from empirical equations; however, it could not simulate reversed flows [11]. Tauveron proposed new empirical correlations to estimate the steady compressor performance in all regions and carried out the dynamic simulations [12].
Most of these models rely on steady compressor characteristics for the entire flow regions, which can only be obtained from experimental data or derived from generic performance curves. Therefore, for preliminary design purpose, it is worthy to use empirical correlations to predict the physical trends of compressor steady characteristics and use them further for transient simulations. Such a 1D model can be regarded as a good compromise between the accuracy and computational load for dynamic simulations in multistage compressors. It can provide more accuracy and more details of the flow field than the lumped method and requires much less computational resources than 3D calculations.
3. Model Description
It is assumed that the flow path meanline is a streamline and that the flow is axisymmetric. The fluid is treated as an ideal gas. The quasi-1D Euler equations are applied in a curvilinear coordinate system. The mass flow, meridional momentum, circumferential momentum and energy equations can be described in vector, nondimensional, conservative form, as follows: The cross-section area is calculated through the simple form .
The source terms represented by in the Euler equations are used to represent the influence of the blade rows and the flowpath on the fluid. They can be splitted into four contributions: the isentropic action of the blades , the influence of friction , and the influence of the flowpath geometry :
3.1. Friction Source Term
In the Euler equations, the entropy production is due to a separate distributed friction force in the momentum equations. This friction force is aligned with the relative velocity but has the opposite direction:
The entropy production can be written as where is the relative velocity vector. The magnitude of the friction force can be rewritten as
Only the tangential component of the friction force can change the energy balance, and the friction source term can be described as following:
3.2. Blade Source Term
This term is implemented to represent the isentropic influence of blade rows. An angular momentum balance over a blade row provides the tangential component of the corresponding (inviscid) blade force :
It is intended that this force is applied in an isentropic way, the viscous contribution being applied separately through the friction source term. Therefore this blade force is perpendicular to the relative velocity:
Once again, only the tangential component of the blade force appears in the energy balance:
3.3. Geometrical Source Term
This term represents the effect of the variation of the cross-section area and the meanline radius. The axial component can be derived as
This geometrical force does not contribute to the energy balance and also does not create entropy.
In the rotating stall and reversed flow regions, the axial velocity may become extreme small, especially around the zero mass flow point. It will generate unphysical magnitudes of both blade and friction forces when they are decomposed in the axial and tangential directions. In this case the blade and the friction forces are combined together and derived from the steady Euler momentum equations applied between the inlet and the outlet of the blade row. These forces and the corresponding shaft work are distributed linearly across the control volumes inside the blade row. The geometrical source term is kept unchanged.
4. Empirical Correlations
The accuracy of the simulations is highly dependent on the source terms in the equations, which are derived from empirical correlations, particularly loss coefficient and deviation angle models. Different empirical correlations have been chosen for all flow conditions, such as normal forward flow, stalled flow, and reversed flow (Table 1).
The separation of boundary layers in compressors has not been fully understood. There are several basic criteria to detect it, such as the static-total pressure coefficient, incidence angle, or the diffusion rate. All these methods are based on experimental observations and their accuracy is highly dependent on the configuration of the compression system. In order to eliminate the uncertainty in numerical simulations, the stall inception and the stall recovery point are obtained from experimental results so that only the loss and deviation models can affect the steady performance. According to the experimental results of Day [2], the stalled point is imposed when the mass flow coefficient is less than 0.45 and the recovery point from stalled flow is set when the mass flow coefficient is more 0.50. The reversed flow is defined when the average axial velocity in the compressor duct becomes negative.
5. The Equation Solver
The computational domain is the real compressor flow path and discretized by a one-dimensional structured mesh. The mesh cells are distributed not only in the nonbladed regions but also inside of blade rows, to make sure that the flow deflection is smooth.
For each control volume, the integral form of the 1D axisymmetric Euler equations is written with the conservative variables located at the cell centers. A quadratic reconstruction scheme with a hybrid limiter is applied to extrapolate the variables at the interfaces. Numerical fluxes are based on Roe’s flux difference. Both explicit and implicit schemes are implemented. The explicit scheme is based on a fourth-order Runge-Kutta algorithm. The implicit scheme is based on an Euler single-step forward in time differentiation. This implicit time-marching method is used for stationary simulations, allowing large local CFL numbers to speed up the calculation.
6. Steady Simulations
This model has been applied to a low-speed multistage axial compressor studied by Day. The test rig consists of four repeating stages with an inlet guide vane and a constant cross-sectional annulus area. The hub/casing radius ratio is 0.75, and the rotation speed is 3000 rpm, giving a maximum Reynolds number of 1.7 × 105, based on the chord [2].
The steady compressor characteristic is simulated using the current model and is compared to the experimental results. The empirical correlations are used to provide loss coefficients and deviation angles for each blade row in all flow regions, such as normal forward flow, stalled flow, and reversed flow. During the numerical simulations, a throttle is gradually turned down in order to obtain the different operating points. When the throttle is set across the instability point, the compressor operates in the stalled region. Meanwhile the reversed flow is obtained by imposing lower static pressure for inlet boundary, stagnation pressure, and temperature for outlet boundary. The comparison of the numerical results and the experimental data are presented in Figure 1.
It can be seen from Figure 1 that the current model underestimates the compressor performance in all flow regions; however the trend is correctly captured. It is found that a discontinuity exists at zero mass flow for both the experimental and numerical results, which is due to the radial circulation [18]. It should be emphasized that all empirical correlations introduced into the model come from the open literature. As far as these correlations can be replaced by a manufacturer’s own correlations, this simulation tool would certainly be able to provide better predictions. At last, it should be pointed out that the purpose of empirical correlations is to provide steady compressor performance that is required for the transient simulations. Once chosen, they do not affect the transient simulation results for a variety of compressor/plenum configuration.
7. Transient Simulations
The time needed for stall cells to be fully developed can be several rotor revolutions and the mass flow undergoes significant changes during this time. The blade forces extracted from empirical correlations are steady and are not adequate for transient simulations. Therefore a first-order time lag equation has been included in the model in order to provide the dynamic blade force [3]. It should be emphasized that only the axial force is modified. The shaft work, and the tangential force are assumed constant
The time constant is approximated to be two rotor revolutions and will be kept constant in all transient simulations. The time lag equation is used only during rotating stall. The experimental observations suggest that the transition between rotating stall and reversed flow does not really occur at zero mass flow [18]. Thus the time lag equation is applied until the flow coefficient reaches the value of −0.09; the steady blade force is used during the rest of reversed flow. The flow range for the application of the time lag equations is shown in Figure 2.
A throttle is imposed at the outlet according to [19] The throttle is further turned down with larger coefficient .
A characteristic method is implemented for the boundary conditions in order to avoid numerical spurious reflections. For positive flow operation the stagnation pressure and temperature are imposed as inlet boundary conditions, and ambient pressure is applied at outlet boundary. In reversed flow operation, it seems unreasonable to expect surrounding air to flow back into plenum through throttle and ambient pressure is always fixed at outlet boundary. Meanwhile, it is likely that the flow can be discharged from the compressor through the inlet boundary. Therefore the ambient pressure is imposed at the inlet when the inlet average velocity becomes negative, otherwise the stagnation pressure and temperature are imposed.
A fourth-order Runge-Kunta algorithm is used for time integration in order to obtain time accurate transient performance. The plenum and throttle are added to form the compression system. A typical mesh is shown in Figure 3.
is used as a reference parameter to determine whether rotating stall or surge may occur in a particular compression system. Different compressor/plenum geometry configurations have been tested in the numerical simulations.
8. Rotating Stall
Although rotating stall is a 2D phenomenon which includes a strong nonaxisymmetric distribution of mass flow and pressure, the present model is able to predict the overall performance when the rotating stall is maturely developed. Since empirical correlations are applied to calculate annular average performance through blade rows, average compressor performances can be predicted by the current model.
The plenum is designed to be 2.41 cubic meters and the compressor rotating speed is fixed at 3000 RPM, and the corresponding parameter is equal to 0.38 [2]. In the numerical simulation, the compressor operates in the normal flow region and it moves into rotating stall as the throttle is gradually turned down. In the experimental test, the compressor operates in a stable manner near the peak of pressure and the throttle is slowly turned down to cause the instability. Figure 4 shows the plenum pressure coefficient versus time (rotor revolutions) for both experimental results and numerical results. It can be seen that the plenum pressure oscillates several cycles before it reaches a stable condition. In the experimental test, the plenum pressure almost keeps constant before the collapse because the compressor starts around the maximum pressure ratio, while the plenum pressure increases in the numerical test since the compressor operates far from the peak point. This is confirmed by Figure 5, which shows the trajectory of the compressor operating point during rotating stall. As the compression system becomes unstable, the operating point settles down on the fully stalled characteristics after a large mass flow excursion. The operating point does not cross the reversed flow neither the recovery point, thus deep surge or classic surge are excluded. This is basically due to a small plenum volume which can be emptied or filled in a short time.
9. Surge
The plenum volume is, respectively, increased to 3.85, 7.61, 12.80 cubic meters in order to obtain parameter to be 0.48, 0.67 and 0.88 at the same rotating speed 3000 RPM. Figure 6 shows the numerical and experimental surge cycles including the steady state curve for the three different plenum volumes. The surge trajectory follows the steady reversed flow curve for a short period. The excursion along the reversed steady curve decreases as the plenum volume is reduced. The smaller the volume is, the faster the flow can be discharged through the throttle and the compressor, and the surge trajectory intersects the reversed curve at lower pressure ratio. For the same reason, the surge trajectory recovers to higher pressure ratios with smaller plenum volumes. Another interesting point is the overshoot of pressure ratio between the surge trajectory and the steady curve around the peak point. It is due to the dynamic deceleration of the flow during the surge, which means that the higher plenum pressure is used to decrease the mass flow in the compressor duct.
In order to observe the flow mechanism in the compression system during the surge cycle, the flow parameters are displayed for the compression system . Figure 7 shows the numerical and experimental static pressure differences between the compressor exit and the plenum inlet. Figure 8 shows the mass flow coefficient at different locations. The static pressure difference between compressor exit and plenum inlet can be seen to accelerate or decelerate the compressor flow during the surge cycle. The deceleration 1 takes place when the compressor stalls (curve AB in Figure 8) and indicates a fast drop of the mass flow. The deceleration 2 takes place when the plenum is repressurized (curve DE in Figure 8), which means that the mass flow decreases gradually and the compressor operating point rises up along the unstalled part of steady-state characteristics. The static pressure at the inlet of the plenum is somewhat higher than at the exit of the compressor, which explains the overshoot in Figure 6. The acceleration in Figure 7 (curve CD in Figure 8) corresponds to the recovery of the compressor from reversed flow to normal flow. The larger the static pressure difference in Figure 7, the faster the flow acceleration and deceleration in Figure 8. The curve BC in Figure 8 shows that the compressor operates in the reversed mode, which can be seen as a steady process for the compressor. Thus the static pressures at compressor exit and at plenum inlet are overlapped during that portion of time.
The static pressure at the compressor exit exhibits some spikes (Figure 7). Compared to experimental results, these spikes are stronger in simulations, due to the empirical correlations which fail to provide the correct pressure ratio for these regions, but the physical trend is similar. This major spike is caused when the compressor recovers from reversed to forward flow (around point C in Figure 8). It is speculated that the mass flow with high pressure and high temperature stored in compressor inlet moves downstream during the transient process from reversed to forward flow.
This can be demonstrated by Figure 9, which shows the total temperature and static pressure coefficients at three different locations during a surge cycle. The total temperature spike at compressor inlet is generated at rotor revolution 24, by the high temperature fluid stored in the plenum and flowing back during the reversed operation. The work input by the compressor in the flow direction increases the total temperature at the compressor inlet compared to the compressor outlet. At the same time, the pressure at the compressor inlet is significantly increased but still much smaller than at the compressor outlet due to the high resistance of the compressor during the reversed flow. This explains the total temperature spike in the inlet duct (point 1). During the reversed operation, the ambient static pressure is imposed at the inlet as shown in the Figure 9. As the mass flow is discharged through the throttle and through the compressor inlet during the reversed process, the temperature and the pressure gradually decrease in the plenum, and so the temperature and the pressure drop at the compressor inlet (from point 1 to point 2).
When the compressor switches from reversed to forward flow, the high temperature fluid stored in the inlet duct is ingested by the compressor, which can lead to even a higher temperature at the compressor exit reflecting the high temperature flow plus the shaft work (from point 2 to point 3). When the plenum is almost emptied, the compressor gets momentum to repressurize the plenum and the static pressure difference between compressor exit and plenum inlet acts as a force that pushes the compressor to unstalled operating point (Figure 7).
A typical deep surge cycle can thus be divided into four portions: (1) deceleration; (2) reversed flow; (3) acceleration; (4) normal flow region. Figure 10 shows the average velocity at compressor exit during deep surge cycle in Day’s experiments. The proportion of time required for each part is roughly 4%, 28%, 3%, and 65%. In the numerical simulations (Figure 8), the time proportion for each part is about 7.8%, 31.1%, 11.3%, and 50% for . Although it is not clear that what was the parameter B in the experiments, the simulation results are quite similar to the experimental results. Obviously, the reversed flow phase is significantly shorter than the positive flow phase since the plenum is discharged through the throttle and the compressor in the reversed flow.
In the experiments the compressor total duct length was either 2.0 or 3.5 m. In the above numerical results, the compressor effective length was imposed to be 2.0 m. The compressor effective length has been extended to 3.0 m, as well as the plenum volume to keep . Figure 11 displays the surge cycle and the mass flow coefficient versus time. The resulting surge cycles have almost the same shape, the parameter playing a predominant role. However, the frequency becomes smaller as the length increases. This can be explained by Figure 12, which shows the experimental frequency, numerical frequency, and Helmholtz frequency versus the plenum volume. It was concluded by Day that the compressor is more akin to a relaxation oscillation than a Helmholtz resonator-type phenomenon in deep surge [2]. It also can be found that the frequency becomes smaller with larger plenum volume at the same compressor geometry.
10. Conclusions
A reliable prediction tool is developed and implemented to simulate the transient performance of poststall behavior in compression systems. The model is based on an adapted version of the 1D Euler equations with additional source terms, which are estimated by empirical correlations. These empirical correlations are only used to provide the steady state performance and can be replaced by the user’s own ones if the accurate steady state performance is needed.
This model can be used for preliminary design purpose as it is only based on the compressor geometry: the blade angles at mid span and the hub and shroud geometry. The model does neither rely on existing experimental performance curves nor on generic ones.
The comparison of numerical steady curve and experimental steady curve shows that the current correlations can capture the main trends of the compressor characteristics with an acceptable accuracy.
For rotating stall, the numerical results can be compared to the experiments for , and indicate that the current quasi-one-dimensional model has the ability to predict the overall performance of the multistage compression system when the rotating stall is fully developed.
For surge, the shape of the surge cycles and the frequency of the surge oscillations in the numerical simulations are tested against experimental results for , , and . The comparisons indicate that the model correctly captures the surge performance. Meanwhile, it is found that the static pressure difference between compressor exit and plenum inlet accelerates or decelerates the flow during the surge cycle, and the high temperature fluid stored in the inlet duct during reversed flow moves downstream when the compressor switches from reversed to forward flow, which can lead to the quite high temperature flow at the compressor exit. Finally, the influences of the compressor effective length on the surge shape and frequency are studied.
Nomenclature
| : | Speed of sound |
| : | Flow section area |
| : | Nondimensional number |
| : | Specific internal energy |
| : | Force term |
| : | Specific enthalpy |
| : | Throttle coefficient |
| : | Effective length of equivalent duct |
| : | Pressure |
| : | Source term |
| : | Mass flow |
| : | Radius |
| : | Entropy |
| SW: | Shaft work |
| : | Temperature |
| : | Rotor velocity at mean radius |
| : | Plenum volume |
| : | Absolute velocity |
| : | Relative velocity |
| Vol: | Blade row volume |
| : | Density |
| : | Constant time scale |
| : | Rotating speed. |
| 0: | Ambient condition |
| 1: | Blade leading edge |
| 2: | Blade trailing edge |
| : | Meridional direction |
| : | Tangential direction |
| : | Plenum |
| : | Compressor. |
| 0: | Total condition. |