Sensitivity Study of Greitzer Model Based on Physical System Parameters of Radial Compressing Units

: Centrifugal compressors are key elements of energy systems and industrial installations including ﬂuid ﬂow. Their operating range is strictly limited by the surge phenomenon. The Greitzer model is a known way of simulating the compressor’s behaviour at the surge. In this paper, the parametric study of different versions of the Greitzer model is conducted. There are several versions of this model that include 4 to 2 equation models. The system behaviour depends on the features of the compressor itself as well as of the plenum. In this paper, all terms connected with the compressor were grouped into the “Co” parameter, while those associated with the plenum were grouped into the “Pl” parameter. The study shows how each component inﬂuences the system stability. The comparison of analytical data with experimental results allowed to draw conclusions regarding the way of choosing the model parameters that provide the best simulation of the real system behaviour. The study shows that the system is well simulated by a model with relatively large values of the L c parameter. The length of the compressor parameter L c = 3.57 m was performing well for the machine with impeller radius 0.33 m. Possible explanations of this ﬁnding are presented and compared to the state of the art. This result may provide possible help in adjusting the model parameters for other machines and designing reliable anti-surge systems based on the Greitzer model suited to energy conversion systems.


Unstable Phenomena in Radial Compressing Units
Centrifugal compressors are widely used in a variety of industries. Their key applications include heavy industry, transport, energy sector, home appliances, and others [1]. In most cases, the stable operation of the compressor is crucial for the undisturbed performance of the entire system including connected network and related machinery. The occurrence of the flow instabilities influences the entire system leading to the drop in the efficiency or, in the more severe cases, to the compressor failure causing significant financial losses. Therefore, the research concerning those phenomena and their avoidance is indispensable.
The operation of the radial compressor becomes unstable when the mass flow rate is decreased below a certain critical value. Emmons et al. [2] was one of the first to identify and analyse the unstable phenomena in the centrifugal compressors. An extensive overview of the research on this topic can be found in [3]. Two main unstable phenomena have been identified: the surge and the rotating stall.

•
The flow within pipes is incompressible, inviscid, and one dimensional • Compression in the plenum is isentropic • The temperature is constant in the entire system • The pressure in the container is uniform Compressing system in Greitzer model [27].
The Greitzer model is a system of four ordinary differential equations, where the first three are resulting from conservation laws of mass and momentum for all three components of the system. The fourth equation is responsible for the response delay of the compressor. The Greitzer model [7] in the non-dimensional form reads (subscripts c, t, p denote the system components: compressor, throttle, and plenum, respectively; notation adopted from [27]): where the following coefficients are applied: • Non-dimensional mass flow rate coefficient: • Non-dimensional pressure rise coefficient: • Non-dimensional time coefficient: • Model parameters: • Compressor steady state pressure rise coefficient Ψ c,ss .
with ρ a denoting gas density, U tip -impeller trailing edge speed, ∆p-pressure rise on a component, D 2 -rotor outlet diameter. f h is the Helmholtz resonator frequency [7]: a a denotes the speed of sound, N s -number of rotations of the rotor needed for the rotating stall to form. There were several values used in the literature including N s = 0 [8], N s = 0.5 [9] and N s = 2 [11]. N s = 0 is, in fact, equivalent to not including the equation in the model. In the case of radial compressors, it is a common practice to omit this equation and limit the applicability of the model to the investigation of the surge phenomenon only. The rotating stall phenomenon in this class of compressors is more complex and this single equation does not provide a close simulation of its behaviour [11,13,26].
It is a widely used approach to additionally simplify the model by removing the equation describing the throttle dynamics. Such a solution is applied in the situation when the length of the throttle is negligible compared to other dimensions of the system. Then the value of G parameter is insignificantly small, which implies almost instant response of the throttle. As the time of response is infinitesimal, the second equation can be removed [15,28,29]. This approach reduces the number of equations to 2 and the set of model parameters to only B and f h . The Greitzer model is not free from limitations. The meaning of model parameters: B, f h as well as the link between them and the geometry of the system, although theoretically defined, is vague for more complex systems. As reported by van Helvoirt et al. [13] the discrepancies between the physical parameters of the system and model parameters are significant and require further study. The most questionable are values of L c parameters that are often reported to be significantly larger than a physical length of the compressor duct [9].
One of the common solutions of this problem is to fit the model parameters to provide the best reproduction of the surge experimental data [9,12,26]. However, the large power systems cannot surge safely even for a short time. Without the surge pressure signal measured at a given machine, the signal fitting is not available. Therefore, without a deep understanding of the method of model parameters evaluation, the Greitzer model becomes ineffective for surge control.
Even though the Greitzer model is known for many years, there is still not much knowledge regarding the mapping between model parameters and physical parameters of the system. In a recent study conducted by Ng et al. [30] the dependency of the model on its parameters was determined with the Taguchi method. This approach provides an insight into the importance of the parameters governing the model. It is still an open question, how the model parameters are connected with the physical features of the centrifugal compressor and the compression system.
The standard non-dimensional form of the Greitzer model describes the properties of the system with the use of implicit parameters such as B or f h . This approach, although beneficial from the point of view of the generality, combines the physical parameters of different components of the compression system (compressor, plenum, and throttle) into a single parameter. For example, the B parameter depends on the compressor rotor tip speed, length of the compressor, and Helmholtz resonator frequency. The last one is clearly a function of plenum geometry. Considering parameters that describe each component separately would be helpful for a better understanding of the role of component parameters. Moreover, it would help to investigate the behaviour of different components in the analyzed system, for example, the same compressor connected to different ducts and throttles. Therefore, this paper presents a sensitivity study of the Greitzer model with respect to the component parameters rather than implicit non-dimensional parameters.

The Aim of the Study
The goal of this paper is to perform the systematical sensitivity study of the Greitzer model allowing investigating the influence of different elements of the compressor system on the model behaviour. The test includes the most popular variants-the two-equation form and the ones taking into consideration the dynamics of the valve and the response delay of the compressor. The investigation is performed with physical parameters of the compression system rather than the model parameters themselves. Therefore, the following parameters were introduced-Co, Pl and Th. Those parameters group the features of each component of the compression system: compressor, plenum, throttle. Together with the N parameter, they formed the set of sensitivity study parameters. Clearly, the physical parameters are implicitly involved in the model parameters and the non-dimensionalization of operation parameters. Therefore, the results are presented in the dimensional form to capture the real influence of the physical parameters of the system on its physical dynamic behaviour. The obtained results are compared with the experimental data.

Investigated System
The sensitivity study was performed for the DP1.12 blower test stand described by Liskiewicz et al. [31]. It is necessary to consider a particular machine in the sensitivity study to provide the compressor operation curve. The best fit model parameters for experimental data collected at this test stand, computed and presented by Grapow in [26], are used for reference.
It is worth underlining that the system physical parameters such as dimensions are fixed for the considered test stand. In the numerical model, they are, however, varied for the sensitivity study. This is possible and meaningful due to the aforementioned vague connection between physical model parameters used by Greitzer. For example, the L c parameter and the physical length of the compressor do not always hold the same value [13].
The test stand cross-section together with its main dimensions is presented in Figure 2. The air enters the system through the inlet pipe (A). Then the flow is accelerated in the Witoszynski nozzle (B). The impeller (C) is located behind the nozzle. From impeller, the flow enters a vane-less diffuser (D) and after that a circular volute (E). The outlet pipe is ended with the throttling valve. Table 1 presents the most important dimensions of the system.  Table 1. The compressor was driven by an asynchronous AC motor (400 V/15 kVA). Although the design point of the blower was at mass flow rateṁ = 0.8 kg/s and pressure ratio PR = 1.12, the blower was operated at lower rotational speed of f rot = 100 Hz, with nominal mass flow rateṁ n = 0.75 kg/s and pressure ratio PR = 1.08, to avoid impeller damage during the investigation of unstable operation. Although this pressure ratio is relatively low, especially for the radial compressors, in the original experiment by Greitzer [7] a compressor with a pressure ratio 1.07 was used. The operational impeller tip speed was equal to U tip = 103 m/s. The impeller had z = 23 blades, which corresponds to the blade passing frequency of f BP = 2.3 kHz during the operation. The volume of the outlet pipe, being a plenum, was equal to V out = 0.0968 m 3 .
The position of the throttling valve placed at the end of the outlet pipe was described by the dimensionless throttle opening area parameter (TOA). TOA = 0% corresponds to a fully closed valve, while TOA = 100% corresponds to a fully open valve. For the sake of consistency of notation, the same parameter will be used to denote the operation point of the compressor.

Sensitivity Study Procedure
The simulations of the unstable phenomena were conducted using different versions of the model. They will be presented in the following subsections.

Two-Equation Model
The first considered version of the model is the one taking into account both simplifications mentioned in Section 1.2. The following system of two ordinary differential equations describing the dynamics of the compressor and the plenum is obtained: where the dependencies Ψ c (Φ c ) and Φ t (Ψ t ) are defined by the compressor performance curve and the throttle resistance curve, respectively. It is a general approach to model throttle response by a quadratic function [28]: with k being chosen in a way that the curve crosses with the compressor operation curve at the considered stationary point of operation. The inverse of this relation provides the Φ t (Ψ p ) relation. The DP1.12 compressor operation curve was elaborated based on experimental measurements presented in [31]. The applied curve fitting procedure was described in detail by Grapow et al. in [26]. It was shown there that the fourth order polynomial has significant advantages over third order polynomials widely used for the estimation of compressor operation curves. The compressor operation curve reads: The dependency of mass flow rateṁ on TOA is given by an experimental relation proposed and described by Liskiewicz in [27]. It reads:ṁ = 0.0943 · TOA 0.5758 (12) The sensitivity analysis was done with respect to two parameters describing the components of the compression system introduced here. The model parameters B and f h depend on a list of system parameters as shown in Formulas (5) and (8), respectively. Some of them, such as impeller trailing edge speed U tip , are well defined and also introduced to the model during non-dimensionalization. Other parameters, such as A c , L c and V p are defining the geometrical model of the compressing system shown in Figure 1. Mapping of those parameters to the real system is to be investigated here. Therefore, the two separate parameters defining the geometrical properties of the simplifies compressing system are defined as follows. The first one is the compressor parameter: According to [7], one should consider the influence of the ratio of A c and L c on the system dynamics instead of each of them separately. The second one is the plenum parameter: The parameter Pl is dependent not only on the geometrical properties of the plenum but also on the speed of sound a a . The reason for that is the consistency with the original model derivation from work of Greitzer [7].
The analysis is performed for several points of operation. The system of ordinary differential equations is solved using ODE45 Matlab routine employing explicit Runge-Kutta formula and the Dormand-Prince pair. The results of the model are a subject of quantitative analysis with the use of two tools: averaging of the signal, to evaluate a mean value of the pressure rise or mass flow rate, and the Fast Fourier Transform (FFT) [32]. In a further step, the frequency spectrum resulting from FFT is analysed to find the highest amplitude and the frequency corresponding to it, both for pressure and mass flow rate.
The initial condition for the simulation was defined as the currently considered point of operation (crossing point of the throttle and compressor curves). To assure that the initial condition is not influencing the result, we performed the FFT analysis only on the stabilized solution. It can be easily shown that the Greitzer model solution stabilizes at the same periodic solution despite the choice of the (feasible) initial condition.

Three-Equation Model-Compressor Response Delay
The second model analysed here takes into consideration the compressor response delay. The model reads: where the dependencies Φ t (Ψ t ) and Ψ c,ss (Φ c ) are defined by Equations (10) and (11), respectively. The sensitivity analysis is performed in a similar manner as in the case of the two-equation model, with respect to Co and Pl parameters, for different values of N parameter, which is involved in η (7).

Three-Equation Model-Valve Modelling
The third model takes into consideration the valve dynamics instead of the compressor response delay. It reads: with the dependencies Φ t (Ψ t ) and Ψ c (Φ c ) defined by Equations (10) and (11), respectively. The sensitivity analysis is performed in a similar manner as in the case of the two-equation model, with respect to Co and Pl parameters, for different values of Th parameter defined as: This parameter determines the dynamic behaviour of the throttle. The parameters A t and L t are involved in G parameter (6).

Comparison with Experiment
The aforementioned test stand was equipped with two dynamic sub-miniature Kulite transducers connected to the Iotech Wavebook 516/E data acquisition system. One of them was placed at the inlet of the system and the other at the volute outlet of the compressor (see Figure 3). The data were collected with sampling frequency equal to 100 kHz. The measurement period was equal to 20.97 s (2 21 samples).
The experimental signals were treated with the FFT algorithm to find the highest amplitude and the frequency corresponding to it. The average value of the pressure signal was not taken into consideration as it was used to determine the approximate compressor operation curve.

Two-Equation Model
In this section, the results obtained with the use of the two-equation model are presented. Figure 4 presents the maps of the dependency of the amplitude of mass flow rate oscillations in the system as a function of Co and Pl parameters for three different valve positions. Figure 5 presents the amplitudes of pressure rise in the system in a similar manner. Figure 6 presents the map for the frequency of oscillations in the system, it is not distinguished here between the mass flow rate and the pressure rise because it can be easily shown that the frequency of oscillations is the same (up to truncation error) for both state variables. The white region in the plot corresponds to the region of negligible oscillations resulting in meaningless frequencies. The Figures 7 and 8 present the maps for mean values of the pressure rise and the mass flow rate, respectively.

Three-Equation Model-Compressor Response Delay
In this section, the results of the model with compressor response delay are presented. From now on, only the pressure rise is considered. There are two reasons for that. First of all, dynamic measurements of pressure rise are significantly easier to perform on the experimental rig than the mass flow rate measurements. Secondly, as can be seen comparing Figures 4 and 5, the nature of the dependency on Co, Pl, and TOA is similar for those state parameters. Figure 9 presents the map of the dependency of the frequency of oscillations on Co and Pl for different values of N. Figure 10 presents a similar map of the amplitude of pressure rise oscillations. The white region of the plot corresponds to the region of negligible oscillations resulting in meaningless frequencies. All results were calculated for TOA = 4.24 (almost fully closed valve).

Three-Equation Model-Valve Modelling
In this section, the results of the model with valve dynamics modelling are presented. Figure 11 presents the map of the dependency of the frequency of oscillations on Co and Pl for different values of Th. The white region of the plot corresponds to the region of negligible oscillations resulting in meaningless frequencies. Figure 12 presents a similar map of the amplitude of pressure rise oscillations. All results were calculated for TOA = 4.24.

Comparison with Experiment
This section presents the results of the comparison of the two-equation model with the experimental data. This model was selected for the comparison with experiment because the number of its parameters coincide with the number of features extracted from the experiment. The model has two parameters: Co and Pl. The pressure oscillations signal recorded on the test stand, after the Fast Fourier Transform, provides two features: the amplitude of the strongest mode of oscillations and the corresponding frequency.
The comparison with the experiment was performed for TOA = 4.24 as it was expected to provide the most reliable data. It was selected because this was the most closed valve position for which the measurements were feasible. In the case of more open valve positions the oscillation modes resulting from surge becomes less dominant with respect to other unstable phenomena, not simulated by the Greitzer model (inlet recirculation or partially rotating stall). If the valve is opened even more, the stable operation of the compressor begins and the model is not producing informative data anymore. The Figure 13 presents the sensitivity maps identical to those presented in Figures 5  and 6 with the horizontal cut at the value obtained empirically. The cut lines for pressure oscillations amplitude and frequency are presented together in the form of a contour plot of the absolute value of the difference between numerical and experimental values. The optimal model parameters Co and Pl are determined by the intersection of those lines.

Two-Equation Model
In Figures 4 and 5 one can see that the amplitude of oscillations depends on the ratio of Co and Pl parameters. There is clearly a border value of this ratio beyond which the system starts to operate in an unstable manner. This ratio corresponds to the critical value of B, where the bifurcation is observed [33]. Please note that the systems with the same compressor, characterized by Co parameter, but connected to plenums with different volumes, will behave differently. The valve position (TOA) significantly influences the positioning of this borderline as well as the severity of the oscillations. The more open the valve, the larger the stability area (the blue region on the map).
In Figure 6 one can see that the frequency of oscillations grows with the growth of the system parameters Co and Pl. The frequency map is not significantly influenced by the valve position. The dependency of the frequency of the oscillations in the system on both parameters on top of the dependency of amplitudes on their ratio assures the uniqueness of the set of parameters resulting in the given oscillatory behaviour.
In Figures 7 and 8 it can be seen that the average value of the pressure rise and the mass flow rate is not equal to the one defined as the operation point in the unstable regime. This means that the occurrence of the instabilities changes the average value of the state variables. Therefore, it can be stated that the oscillations are not symmetric with respect to the operation point. Monitoring the average values of state variables can be considered to be a way to control the stability of the machine. On the other hand, this means that measuring the state variables in the unstable regime and averaging them to determine the operating curve of the compressor in that area, although widely done, may lead to significant errors.

Three-Equation Model-Compressor Response Delay
In Figure 10 it can be seen that the growing N parameters enlarge the stability zone. It deflects the stability line to the larger Co values. The amplitude is no longer dependent purely on the ratio of Co and Pl parameters. In Figure 9 one can notice that the frequency of the oscillations does not depend significantly on the N value. N = 0 corresponds to the two-equation model, to avoid numerical errors the simulation was done for a small non-zero value N = 0.01 and used for reference.

Three-Equation Model-Valve Modelling
In Figure 12 it can be seen that for the decreasing values of the Th parameter, the stability zone in the map is shrinking. The borderline between stable and unstable operation is deflected towards lower values of Co parameter. In Figure 11 one can notice that the frequency of the oscillations does not depend significantly on the Th parameter. Th = 1 corresponds to the two-equation model, to avoid numerical errors, the simulation was done for Th = 0.5 which return very similar results. Those results were used for reference.

Comparison with Experiment
In Figure 13 one can see that the parameter values resulting in the same amplitude and frequency of pressure rise oscillation as the one obtained during the experiment are Co = 4.2 × 10 −3 and Pl = 1.15 × 10 6 .
According to Equation (14) and assuming the speed of sound at the operating condition to be 343 m/s and the volume of the outlet pipe from Table 1, we can calculate the theoretical value of the Pl parameter. This is equal to Pl th = 1.21 × 10 6 . The numerically obtained value is very close to the theoretical one. The discrepancies can be caused by the shape of the considered plenum volume. In our experiment we used a pipe with a length significantly larger than the diameter, not a container as assumed by Greitzer [7].
According to Equation (13) the Co parameter is just a ratio of the area of the compressor to its length. The meaning of those parameters is vague. If we assume that the area of the compressor parameter in the model is the area of the compressor outlet, which can be calculated using the data from Table 1, then A c = 0.015 m 2 . In such a case, to achieve Co = 4.2 × 10 −3 , the length of the compressor must be equal to L c = 3.57 m. Looking at Figure 2 and Table 1 one can notice that this value is significantly larger than the length of the duct through the compressor. If we consider any other cross-section of the duct to compute the area, the results are not significantly different. According to [7], the A c L c ratio (here denoted as Co) should be computed as an integral along the duct. This, however, will not result in the desired change of the physical value of Co parameter due to insignificant (less than an order of magnitude) variations of the duct area. One of the other reasons for this discrepancy could be a relatively low pressure ratio of the considered machine preventing the occurrence of more severe forms of surge. On the other hand, similar large values of L c were reported in the literature before for even smaller radial machines, Table 2 presents an overview of the values of L c parameter compared with the radius of the impeller for several Greitzer model applications presented in the literature. It is clearly visible that L c values are significantly larger (by at least an order of magnitude) than radii of the impellers, which is not physical. Apart from that, Van Helvoirt et al. [13] explicitly reports the necessity of increasing L c value from 0.197 m to 1.8 m and 3 m for the meaningful simulation. Let us investigate if the additional equations with the feasible selection of parameters could allow more realistic selection of L c parameter. In the case of the throttle dynamics equation the selection of small values of Th is causing the deflection of the stability limit line to the left as presented in Figure 12. This would cause the yellow line in Figure 13 to deflects in similar manner if this model would be applied. The smaller value of Co and therefore even larger L c would be obtained. Note that the values of Th can only be positive and smaller than 1. Therefore, there is no feasible value of Th parameters that can minimize the discrepancy between L c and the real size of the compressor. In the case of the compressor delay equation the selection of N larger than one results in deflecting the stability limit line to the right as presented in Figure 10. This seems promising. Further investigation, however, shows that the frequency plots of the two equation and the compressor delay models (Figures 6  and 9, respectively) do not differ significantly. The main difference is that the frequency values for the compressor delay model are slightly higher. Therefore, the blue line in Figure 13 would be lying closer to the axes of the plot. Adding to that the deflection of the yellow line to the right would result in slightly better fit in terms of L c (higher Co, meaning lower L c ). On the other hand, a good fit of the Pl and the plenum volume would be lost. Therefore, adding the compressor delay equation is not a cure either.
One can conclude that the very high value of the L c parameter has no connection with any actual size of the compressor. This leaves a question of whether the L c should be purely connected with the geometrical parameter of the machine as it is often claimed [34]. In the radial machines, the radial velocity component of the gas is significantly larger than in the case of axial compressors considered originally by Greitzer [7]. Additionally, the gas is present in the rotor for a significant amount of time (larger than in case of axial compressors), even for several revolutions. This violates the assumptions made by Greitzer. It also suggests that the A c and L c parameters might need to be approximated differently for axial and radial impellers. This topic should be investigated further to allow efficient modelling of the compressing system with the use of the Greitzer model without access to the experimental data.

Summary
This paper presents the sensitivity study of the Greitzer model in its two-equation form as well as the versions taking into consideration the delay of the compressor response and the dynamics of the throttle. It presents the influence of several physical parameters of the system on its behaviour. Three parameters Co, Pl, and Th were introduced to represent features connected with each component of the compression system (compressor, plenum, and throttle). Together with the N parameter, they formed the sensitivity study parameter set. Such an approach, considering the physical parameters connected with particular system components, instead of standard implicit non-dimensional parameters, provides a better understanding of the role of each component. It also allows us to consider the behaviour of different components in the system, for example, the same compressor connected to different ducts or throttles.
The analysis was done on the example of the DP1.12 blower. The results obtained with the numerical model were compared with the experimental data collected for the considered machine. There was an excellent match between the experiment and model in case of the plenum volume. The physical parameter describing the properties of the compressor had a value significantly different than expected. Namely the length of compressor parameter L c = 3.57 m was found suitable for the machine with impeller radius 0.33 m. The analysis of other studies presented in the literature has shown that such a situation is not unusual.
It was shown that incorporating considered additional equations will not significantly decrease the discrepancy between theoretical and determined values of Co parameter without worsening the fits of other parameters. Additionally, expanding the model by additional equations with uncertain coefficient values seems inexpedient. In the future, the presented framework of fitting the parameters and the conclusions concerning the applicability of additional equations can be used for an analysis of a variety of compressors to draw more general conclusions.