Inverse Design of Centrifugal Compressor Stages Using a Meanline Approach

This article discusses an approach for determining mean-line geometric parameters of centrifugal compressor stages givenspeciﬁedperformancerequirements.Thisiscommonly knownastheinversedesignapproach.Theoppositeprocess, thatofcalculatingperformanceparametersbasedongeom-etryinputisusuallycalledanalysis,ordirectcalculation.An algorithmandcomputercodeimplementingtheinverseap-proachisdescribed.Asanalternativetocommerciallyavail-ableinversedesigncodes,thistoolisintendedforexclusive OEMuseandcallsatrusteddatabaseoflossmodelsforindi-vidualstagecomponents,suchasimpellers,guidevanes,dif-fusers,etc.Thealgorithmextendsapplicabilityoftheinverse designcodebyensuringenergyconservationforanywork-ingmedium,likeimperfectgases.Theconceptoflosscoefﬁ-cientforrotatingimpellersisintroducedforimprovedloss modelling.Thegoverningconservationequationsforeach componentofastagearepresented,andthendescribedin termsofaniterativeprocedurewhichcalculatestherequired one-dimensionalgeometry.Agraphicaluserinterfacewhich facilitatesuserinputandpresentationofresultsisdiscussed brieﬂy.Theobject-orientednatureofthecodeishighlighted asaplatformwhicheasilyprovidesformaintainabilityand futureextensions.

lyzing geometry rather than doing inverse design.Furthermore, many of these codes have no generalized loss models, but have instead past compressor performance data tabulated to indicate performance of the machine under study.This technique works well for cases where the new design is similar to an existing one, but poorly when the new design differs significantly.
These codes are used for design in an iterative manner, that is, the designer inputs geometry and sees if the code predicts the desired performance.A final design is thus arrived at when the desired performance is arrived at.A more direct route to the final design was thus desired.
Several commercial codes are available which solve the design problem directly, such as Concepts NREC's Compal program (ConceptsNREC, 2002) or PCA Engineers' OONA program (PCA, 2002).However, these codes use a generalized set of loss models which are not specifically applicable to this OEM's test experience.It was desired to incorporate the basic idea of these codes with the specific loss model capability already present at our facility.
It was decided not to pursue investment in a commercial code for other reasons as well: a) the new code was to serve as a platform for future extensions into 2-D and 3-D analysis, b) extensive modification of the code was seen as a natural part of its evolution, and c) a code with a GUI written with a multiplatform library was required.To meet these objectives, a highly object-oriented code was written with a flexible computational engine.This would allow drop-in-place changes such as extensions to equations of state, loss models, and other geometric components (e.g.sidestreams).It was felt, therefore, that a code which the OEM has detailed knowledge of, and control over, would best meet these objectives.

ALGORITHM
The following calculation procedure is based on well known 1-D equations for centrifugal compressors.Indeed, most of them can be found in textbooks such as Dixon (1998) or Japikse (1996).Of interest here is the manner in which these equations are solved to provide a unified design approach for the entire stage.The user begins by selecting the inlet total pressure, P t00 , and total temperature, T t00 .Thus any total property at the inlet is known using the equation of state.In particular we are interested in the inlet total enthalpy: h t00 = f eos (P t00 , T t00 ) [1] Also selected is the desired total pressure at the stage outlet, P t80 .By imagining an isentropic process from the inlet to the outlet, we can write the following set of relations: which can be summarized as: This process is illustrated in Figure 1.The program initially guesses a total to total stage efficiency which is used to find the outlet total enthalpy, h t80 .The efficiency used corresponds to the isentropic ideal process rather than the polytropic: This value will be used to find an important parameter for continuing the calculations, the work input to the stage.We will discuss this point momentarily, but for now consider that the computation proceeds by calculating each individual component (e.g., impeller, diffuser, volute, etc.) individually until the outlet state of the final component is known.This state will usually correspond to an outlet total pressure different than that given by the user, and hence, will correspond to a stage with a different total to total efficiency than that which was

FIGURE 2
Overall flowchart of algorithm.guessed above.Thus the program guesses a new stage efficiency and the procedure is repeated until the guessed value equals the calculated value.This step constitutes the outermost iteration loop in the calculation algorithm (see the flowchart in Figure 2).
It is important to note that the individual component calculations are done using external, independent loss models.Within each component calculation, an iteration to find the correct loss coefficient is done.These loss coefficients then determine indirectly the overall stage efficiency.This is the fundamental reason why finding the overall stage efficiency is an iterative process.
The next step in the algorithm is to compute the work input for the stage: Now, the work input coefficient, is given by the user, so the calculation of U 20 can be made in a straightforward manner.Also, either flow coefficient or impeller rotational speed is known from user input.If rotational speed is given, impeller diameter is found from: If the flow coefficient is given, then impeller diameter is found from: At this point all overall parameters required for continuing the calculation have been determined.The procedure continues by calculating the individual components which the user has identified as pertaining to the stage.For the description here, each of these components will be treated in order, beginning with the inlet and ending with a return channel (or return channel vane extension) for an interior stage of a multistage machine, or ending with a volute and discharge cone for a single stage machine or final stage of a multistage machine.Figures 3-4 contain crossectional views of a typical stage.

Inlet Calculation
The design objective for the inlet is finding the shroud diameter at exit (or inlet to the impeller).This is done by minimizing the relative velocity (at the shroud) going into the impeller by varying the shroud diameter (see Figure 5 for a general velocity triangle).The procedure used is to simply guess the shroud diameter, do the calculation, and then determine if the resulting shroud relative velocity corresponds to a minimum.Thus, with the guessed D 08shd a guess is made for the efficiency (static to static).Note that station 07 here is identical to station 00.We use 07 to be consistent with traditional practice.For a vaneless inlet section, these Equations ( 12)-( 13) are replaced by the free vortex relationship and the familiar pythagorean relationship from the velocity triangle: The annular area through which the flow passes is found for later use: Now, if the flow is diffusing, the static-to-static efficiency can be used to find h 08 If the flow is expanding, the equation becomes As before, the outlet pressure is found from the inlet state and the outlet isentropic enthalpy via the equation of state: The outlet temperature and density follows also from the equation of state: Now, a calculated mass flowrate can found, and compared to the known value, ṁ08 .The outlet enthalpy, h 08 is guessed until ṁcalc 08 = ṁ08 .The outlet total pressure can now be found from the static conditions and the total enthalpy: This enables the calculation of a loss coefficient, A loss coefficient can also be found independently through the use of a loss model, The two values of loss coefficient will be equal if the correct value of η 0708 ss was chosen above.This gives rise to another iteration loop to establish the loss coefficient.Now the relative velocity at the exit shroud needs to be found.The impeller rotational speed is known at this point, so the speed of the inlet shroud is found as Recognizing that location 08 is effectively the same as 10, we already know the meridional and tangential velocities (C m10 , C u10 ), thus completing the velocity triangle.The shroud relative velocity is found from: where we additionally assume that the flow angle (β shd 10 ) is the same as the guide vane blade angle preceding it (β bl 08 ).At this point W shd 10 is checked to see if it indeed represents a minimum.If it does not, then a new value of D shd 08 is guessed, and the iteration proceeds.If it does, the inlet calculation is concluded.

Impeller Calculation
The procedure begins with a guess for the loss coefficient, (LC guess ) which constitutes the start of an outer iteration loop.A user input diffusion factor enables the relative outlet velocity to be found as From Equation ( 6), the work input to the impeller can be found.Then, using the Euler turbomachinery equation, C u20 can be found.The following formulas solve the remainder of the outlet velocity triangle: Since the outlet conditions are now known, it remains to find the blade angle through a slip model.The program can adopt a variety of slip models so, in general, For instance, one such formula is the Stanitz model (Stanitz, 1952), Note that here the blade angle rather than the flow angle is needed, and this value has not been computed yet.
To use this model, a guessed blade angle would be needed, a reasonable value being the flow angle from Equation ( 35).
With slip velocity known, the blade angle is determined by finding ideal tangential velocities, as follows: [38] [39] Note that if using the Stanitz slip model above, this would be the new blade angle to be used in the iteration loop.Convergence in this case is achieved very quickly.
With blade angle known (along with user inputs for number of blades and blade thickness), the outlet circumference due to blade blockage can be determined.This value can be used to find the circumference which permits flow passage, and a blockage coefficient: Now, for the impeller only, we are defining the loss coefficient as follows: Using this definition, and remembering that we have guessed a loss coefficient as part of the iterative scheme, we can solve this equation for the total to total efficiency.Then, the isentropic total enthalpy at outlet can be found by solving and the total pressure is found from the equation of state as Since two total properties are known, others follow directly from the equation of state: The outlet static conditions are found as follows: [50] This now permits the calculation of an important design parameter, the outlet width and area: At exit, the flow proceeds into a space whose width is greater than B 20 since there is no blade blockage and because a small gap is present between the impeller and diaphragm walls.The user inputs the relative opening that this gap represents, thus enabling the width at station 21 to be found.Other geometric parameters follow: An isentropic diffusion process is assumed to occur between 20 and 21, thus the total state remains the same: Also, from Equation ( 56) it follows that the tangential velocity component is constant: Thus by solving the following four equations in four unknowns (h 21 , ρ 21 , C 21 , C m21 ), we establish the outlet static state and velocity: Other quantities of interest are found directly from: At this point an impeller loss model is invoked for a new loss coefficient in the same manner as Equation ( 25).This value is compared with the guessed value above and iteration continues until they are equal.The loss model used here was developed specifically for the OEM's impellers.However a number of general models have been introduced which attempt to account for the complex nature of impeller flow.Japikse, for instance, uses a two-zone model where the flow is broken up into a nearly isentropic primary region and a secondary region near the suction surface dominated by wall shear effects (Japikse, 1985).The use of such models in combination with our specific test experience is certainly a subject for future research.

Vaneless Diffuser Calculation
For this description we assume the impeller is followed by a vaneless diffuser, or a vaneless space preceding a vaned diffuser.In either case, the calculation procedure is the same.
In general, the design objective here is to find the diameter and width at the outlet.This is done, as for the impeller, by using a known diffusion factor supplied by the user.However, it is impossible to close the equations without knowing either the width or diameter itself.So, the user is given the choice of setting, in addition to the diffusion factor, either of these parameters.The calculation then solves for the other parameter, the one left as unknown.The description that follows will illustrate the procedure for both cases.
The inlet to the diffuser is taken as the impeller outlet, station 21 for which all geometry, velocities, and thermodynamic quantities are known.As with all our component calculations, a loss coefficient is first guessed.Thus, the following set of equations can be solved directly: [81] To close the iteration, a calculated value of loss coefficient can be found from an external loss model: This value is then compared to the guessed loss coefficient and a new value selected accordingly.

Vaned Diffuser Calculation
The vaned diffuser (either conventional or LSD) begins at station 30 and ends at station 40.Here there are three geometric unknowns of interest, the outlet diameter, width, and angle.In order to close the equations, two need to be selected in addition to the diffusion factor.Thus the calculation essentially solves for either diameter, width, or flow angle depending on what the user selects as known.
As before, the user supplies a desired diffusion factor.A loss coefficient is guessed initially and calculation proceeds in exactly the same manner as shown in Equations ( 73)-( 80), with the proper change of subscripts to reflect the new location.Then, a set of four equations in four unknowns (A 40 , C m40 , C u40 , and one of D 40 , B 40 , or α 40 ) is solved, three of which are the same as Equations ( 81)-( 83).The fourth equation is simply the definition of the flow angle: The loss coefficient is found by using a loss model for the type of vaned diffuser in question.This is represented generally in the same way as Equation ( 85).It should be noted that, as for the impeller, many general loss models exist for diffusers, such as Ribi and Dalbert (2000) and Japikse and Baines (2000).
Note that the vaneless diffuser calculation is used for any vaneless space, including the return bend.Likewise, the vaned diffuser is used to represent all vaned calculations, including the return channel.For this description it is useful to imagine a second vaneless space after the vaned diffuser representing stations 40-50.For a last stage, a volute would begin at station 50 and extend to station 70.

FIGURE 6
Volute and discharge cone configuration.

Volute Calculation
The design objective for the volute is the determination of its outlet width and height.The height is represented in terms of minimum and maximum radii.As shown in Figures 4 and 6, these parameters give the basic size of the volute.In order to close the calculation, a loading coefficient is supplied by the user along with an aspect ratio.As usual, a guess for the loss coefficient is made by the program which initializes the outer iteration of the calculation.A guess for the volute radius is also made, constituting an inner iteration loop.The volute is designated by the user as external (fixed inner radius, increasing outer radius) or internal (fixed outer radius, decreasing inner radius).
The following formulae are by now routine, and establish the total outlet state: ṁ70 = ṁ50 [87] Now, since the loading coefficient expresses a deviation from free vortex flow, by knowing this parameter, the tangential velocity at the exit can be found.
Furthermore, since the meridional component of the velocity at 70 is very small, it is commonly assumed that the overall velocity is equal to the tangential velocity: The following relationships are straightforward and establish the outlet static state.For an external volute the outer radius can now be found: Similarly for an internal volute the inner radius is found as The mean radius is now found: Once again, the loss coefficient is calculated from a loss model and compared to the guessed value: As before, this is a specific, empirically based model, although general models, such as Weber and Koronowski, 1986 are available.

Discharge Cone Calculation
The primary purpose of the discharge cone is to connect the volute exit with the customer's downstream process piping.Generally speaking, the diameter of this piping would be known.The calculation then would find the corresponding velocities and thermodynamic properties.However, in case the diameter is unknown, the user may choose to enter instead a known static pressure, or a known velocity, in which case the diameter is computed by the program.
The calculation of the total properties is the same as that described in Equations ( 87)-( 90), so it will not be repeated here.Note that P t80 is found by the equivalent of Equation 89.Next, the following set of four equations in four unknowns is solved: At this point, it is presumed that the stage outlet has been reached.Notice that a value of P t80 has been calculated independently of the user input value.Because of this, the value of total to total stage efficiency will be different than what was guessed in the beginning as part of the overall stage iteration loop.Therefore, a new total to total efficiency is found based on the actual, calculated conditions at the discharge cone outlet, and the iteration is repeated.

Curvature Calculations
After each component is calculated, the user may invoke a calculation for the meridional velocity profile due to curvature at the hub and shroud.This is based on a discretized form of the Euler-n equation where n is the direction normal to the flow, and θ is the angle of inclination of the section with respect to the radial.This calculation is done after the fact, that is its results have no influence on any downstream component calculation.It simply provides the user with information as to the effects of curvature at a particular section.

Non-Radial Diffuser Calculation
The code provides for the calculation of non-radially oriented diffusers, return bends, and return channels.The radial orientation in this case is only at the inlet and outlet.An additional unknown, the inclination angle, is present in any such calculation although, except for finding geometry at such a section, the calculations remain essentially the same.

Iterative Scheme
Most of the iteration schemes in the foregoing description involve a simple update to the guessed value with the new, calculated value.For cases where convergence is more difficult (such as finding the outlet radius of the volute), iteration proceeds by a) establishing the low and high limit of the solution domain, b) breaking up the domain into an equal number of increments, c) guessing in increments from the initial value to the high limit and then to the low limit, and d) when the solution is bounded, using the bisection method to converge.This scheme, although simpler and slower than other schemes, is extremely reliable.

Equation of State
The reader will notice the general form in which calculations involving the equation of state have been presented.To the calculation algorithm, the equation of state simply represents an external function by which knowledge of a thermodynamic state can be used to find any other thermodynamic property.Thus, algorithmically, it is not important which equation of state is being used.Any equation of state can be written and added without changing the code that does the calculations.Indeed, external libraries of functions with many equations of state can be purchased and linked to the program with negligible modification to the base code.This raises the issue of software design, to which we now turn.

PROGRAM ARCHITECTURE
The algorithm illustrated above was written using the C++ programming language and adheres strongly to the principles of object oriented software design.Thus separate classes are created for each stage component (e.g., impeller, diffuser, volute) and the stage itself, all of these in turn derived from a generic component base class.Likewise, classes are created for gases, equations of state, loss models, iterators, and user interfaces.The idea was to have as much "drop-in-place" capability as possible.So, for instance, if an equation of state should change, the programmer simply has to program a new class for the particular equation of state he is dealing with.No other code modifications would need to be added.
The main user interface for the code uses a GUI library provided by Qt (Trolltech, 2002).This was used to create a menu bar along with several dialog boxes for convenient user input.One of the advantages of this particular library is that it is multiplatform, i.e., it can be compiled under both the UNIX and Windows operating systems.Several support classes were added to the user interface, particularly for allowing unit conversions and manipulating the communication with the calculation engine.It should be noted that, consistent with good software practice, the GUI and calculation engine are completely separate.The GUI only has pointers to the calculation components and the calculation only has a pointer to the GUI.This permits the use of any user interface, such as a command line interface.Any other interface could also be dropped in after the fact.
An emphasis on future maintenance and extensibility was especially required since this code will be used as a front end for 2-D and 3-D flow solvers.Thus, the final product is envisioned as a complete aero-design environment for centrifugal compressors, and many future modifications are anticipated.

OTHER MODES OF OPERATION
The fundamental purpose of the algorithm presented here is to perform a design calculation, that is to find geometry given performance.However, the code can also be run in analysis mode where performance parameters are found given geometry.
A mixed mode of calculation, where a combination of some geometry and some performance parameters are given, resulting in a mixture of geometry and performance, is also possible.This mode is particularly useful for designers who are trying to modify an existing machine and have access to some of the geometric information, and some of the performance conditions.The output is always the full set of geometry, velocities, and thermodynamic properties.
The following describes how each component offers the user the choice of design or performance input parameters.To begin, for the overall stage the user may input either the power or outlet total pressure.The user may also enter both, but if this is done it is equivalent to setting the stage efficiency, so efficiency will not be calculated independently through loss models.
For the inlet (station 00), the user may input, in addition to the total conditions, either the static pressure, crossectional flow area, or velocity (meridional).The same is true for the outlet when it is a discharge cone.
For the inlet guide vane, the user may choose either the angle of the vane, or the preswirl coefficient.Of course the user may also choose a vaneless inlet.
For the impeller the user chooses either the diameter or the work input coefficient, either the diffusion factor or outlet width, and either the mass flowrate or the outlet blade angle.Of the latter, the user may choose both, but doing so is equivalent to setting the slip, so it will not be computed by a slip model.
For a vaneless diffuser, as discussed, the user may choose any two of outlet diameter, width, or diffusion factor.For a vaned diffuser any three of diameter, width, diffusion factor, or outlet angle are chosen.
The basic conservation equations shown above, naturally do not change when the code is run in its various modes.The iterative algorithm does change, however, and must take into account the great many combinations of inputs that a user may wish to run.The details of this will not be discussed here since it is a separate topic.It will be sufficient to note that the various modes of operation were deemed necessary given the diverse needs of the users at the OEM facility.

NOMENCLATURE
FIGURE 4Stage configuration for single stage.

P
70 = f eos (h 70 , P t70 , T t70 ) [94] T 70 = f eos (h 70 , P 70 ) [95] ρ 70 = f eos (T 70 , P 70 ) A 70 and AspectRatio are known at this point, the following two equations are solved for R 70 and B 70 : A 70 = R 70 B 70 [ 80 is given the set of unknowns is (ρ 80 , h 80 , P 80 , C 80 ), if C 80 is given the set of unknowns is (ρ 80 , h 80 , P 80 , D 80 ), and if P 80 is given the set of unknowns is (ρ 80 , h 80 , C 80 , D 80 ).As usual, a loss coefficient calculated from a model is used to compare with the guessed value.Other thermodynamic properties of interest can also be found in a straightforward manner using the equation of state.