Optimization of a Centrifugal Compressor Using the Design of Experiment Technique

Centrifugal compressor performance is affected by many parameters, optimization of which can lead to superior designs. Recognizing the most important parameters affecting performance helps to reduce the optimization process cost. Of the compressor components, the impeller plays the most important role in compressor performance, hence the design parameters affecting this component were considered. A turbocharger centrifugal compressor with vaneless diffuser was studied and the parameters investigated included meridional geometry, rotor blade angle distribution and start location of the main blades and splitters. The diffuser shape was captured as part of the meridional geometry. Applying a novel approach to the problem, full factorial analysis was used to investigate the most effective parameters. The Response Surface Method was then implemented to construct the surrogate models and to recognize the best points over a design space created as based on the Box-Behnken methodology. The results highlighted the factors that affected impeller performance the most. Using the Design of Experiment technique, the model which optimized both efficiency and pressure ratio simultaneously delivered a design with 3% and 11% improvement in each respectively in comparison to the initial impeller at the design point. Importantly, this was not at the expense of sacrificing range, of critical concern in compressor design.


Introduction
Design and optimization of rotating machines by means of the numerical optimization methods has become an established technique in all industries for many decades.
In centrifugal compressors, the impeller plays a most important role in machine performance and this has attracted significant interest of researchers' and developers over several decades. Not only will a good impeller design bring significant performance benefits itself, but a good impeller will deliver a more uniform flow to the diffuser hence improving the performance of the whole compressor.
In the preliminary design stage, inlet and outlet angles and the main dimensions of the impeller are obtained. This design will be followed by detail design where the complete gas path and the geometry are specified which includes hub and shroud contours and the blade curvature for the vaned components [1]. During detail design some adjustments are repeatedly necessary as the acceptable performance in a relatively wide design range is constrained by some limitations [2]. Thermodynamics and fluid mechanic principles, mechanical, vibration, and the manufacturing constraints and the cost dictate these limitations. Design and optimization of rotating parts by means of the numerical optimization methods has become an established technique for recent decades. Numerical optimization techniques help compressor designers obtain optimum efficiency or maximum pressure ratio [3] or widener stable operating range [4].
Evaluating the effect of different aspects in the optimization process simultaneously imposes heavy computational cost, hence some techniques are developed to reduce the problem exploring space. The Design of Experiment (DoE) technique is one of this methods. The DoE technique helps the designer limit the calculations needed to find the effects of input variables on the output variables. In this technique the shortcoming of inaccuracy of using an approximation function is compensated by reducing the number of calculations. A combination of mathematical and statistical methods are used for constructing the design space [5].
Many researchers have targeted the optimization of centrifugal compressors. Many parameters are considered in this regard. For example optimization techniques were implemented to optimize the blade profile [6], meridional profiles [7,8], blade sweep angle [9], the splitter [10], the vaneless diffuser [11] and there are many others.
In the metamodeling technique to build a surrogate model, DoE can be employed to generate appropriate training samples [12]. Barsi et al. optimized the geometry of a return channel for a centrifugal compressor through the artificial neural network (ANN) method which is trained using a database of 3D CFD solutions created with a DoE technique [13]. Kim et al. also optimized the meridional geometry of an impeller. They used Latin-hypercube sampling of DoE to generate design points for implementing the Neural Network (NN) method [14].
Siemann and Seume designed a compressor diffuser using DoE to see the effects of four parameters and their interactions on the blade efficiency [15].
Bonaiuti et al. presented a parametric study for an optimization problem using DoE. They identified the most influential geometrical factors affecting impeller performance and optimized three different impellers [16]. Guo et al. optimized a centrifugal compressor impeller using DoE and response surface models to increase pressure ratio, isentropic efficiency and to reduce work input of a compressor [17].
Tang et al. introduced a response surface method for constructing a surrogate model to parametrize the design space for aerodynamic optimization of a transonic fan [18]. Wang et al. considered the β angle distribution to optimize a centrifugal compressor impeller. They implemented the Kriging interpolation to construct the approximation function [19].
In this article three main characteristics of an impeller are considered for investigation of their effects on performance. The meridional contour geometries, the angle distribution of the blades and the leading edge of main blades and splitters are parametrized in this regard. Parameterization of the meridional geometry also defines the vaneless diffuser including the level of pinch, vital ensure capture of the full compressor geometry. The Full Factorial method is used in each step to explore the most influential factors. The models are created and solved by implementing a three dimensional viscous solver and the CFD outputs are used to obtain each model performance including the pressure ratio and isentropic efficiency. The Box-Behnken method is used to construct the design space for the selected parameters. Second-order polynomial regression is used to construct the response surface for optimizing pressure ratio and efficiency simultaneously.

Bézier Curve
Bézier curves are selected for describing the geometry of curves. A Bézier curve of degree n is a polynomial curve which is expressed in the following form, in which, the c i are the control points and B n i is ith Bernstein polynomial of degree n, by following definition Bézier curves allow control of the curves systematically and the derived curve is smooth with continuous derivatives [2]. The lines between consecutive control points are named the control polygon. The Bézier curve passes through only the first and the last control points and tangent to the control polygon at both ends. The most influential points are the end points and the middle control points have lower weight in creating the curve and can be used for modification the curve. Figure 1 shows the effect of changing a control point on a Bézier curve.

Design Space
To investigate the effect of each factor on the objective function, the full factorial design technique can be used. A 2 k full factorial design space is built by k parameters each at two levels results in 2 k design points. This type of design is helpful in early stages of design to recognize the most effective parameters in the results.
In the case with A, B, C and D as the factors, each at two levels, 2 4 combinations are generated. Four of them are the main effects shown by A, B, C and D estimate the main effects of the related factors on the the result. In six of them, the two factor interactions considered are shown by AB, AC, AD, BC, BD and DC. The three parameters interactions are ABC, ABD, ACD and BCD. The four factor interaction for this case is ABCD. The main effect of A, is the average of the results, in which the combinations with A at the high level are positive ones and the combinations with A at the low level are negative ones. For each parameter two levels are shown by + or − respectively. If the average response at the combinations in which X at the high level is noted by y x + and the average response with X at the low level is noted by y x − , The effect of X combination on the objective function [5] will be, To recognize the significant factors, a normal or half-normal probability plot analysis is the simplest tool. Normal probability plots show that whether sample data conform a normal distribution or not [5].
For constructing a normal plot, the effects are first ranked from smallest to largest, E 1 < E 2 < ... < E n , Then the ordered E i are plotted versus: Then a straight line interpolated between the 25th and 75th percentile points. The factors that have negligible effect on the objective function have normally distributed Es and they fall along a straight line, whereas most effective factors are far from the straight line and hence normal distribution [5].

Box-Behnken Technique
Different types of designs are being used for creating the design space. The Full Factorial Technique is a useful technique for exploring all of the design space. With n parameters, each at k level, the full factorial creates k n design points. Two-level factorial techniques are very efficient but they result in first-order models and cannot detect curvature. At Three-level, the 3 n factorial, is not efficient and generates huge design space points [20].
Box-Behnken is a technique which makes an efficient design space and it requires only the fraction of full factorial points to estimate the second-order effects of the response surface. For A, B and C factors each at three levels {−1, 0, 1}, Box-Behnken design offers the points as shown in Figure 2 which can be compared with 27 design points obtained from the full factorial technique.

Response Surface Model
After building the design space for each combination, the response (objective function) could be derived. The response surface is constructed as a surrogate model based on the most important effects. For instance if X,Y and XZ combinations show the largest effects on the objective function, the estimated response surface [21] is defined as follows, in which β 0 is the average of the response for all combinations, β x is the E X /2 and so on. In this relation the X,Y and XZ factors are varied in [−1, +1] intervals which correspond to variations between low and high levels.
To find the relationship between the objective function and factor values, the following second-order polynomial model is used, in which {x 1 , ..., x k } are the factors and β i s are the coefficients and is the approximation error.
The low order polynomial models may not be a reasonable approximation for entire design space; however they work quite well over a relatively small part of the domain. The method of least squares is used to estimate β coefficients in Equation (6). This method chooses β i s coefficients so that the sum of the squares of the errors, L, is minimized.

Creating the Design Space
The parameters for the first three stages of the design space created are shown in Table 1. The 2 k full factorial design technique is implemented in these stages for creating the design space. After performing stage I to III the most effective parameters are recognized for the design space created in stage IV.  Figure 3. A hub curve control polygon was constructed between two fixed points H 1 and H 2 and two control parameters M 1 and M 2 . Fixed points and control points for constructing a Bézier curve for shroud profile are S 1 , S 2 and M 3 , M 4 respectively. The diffuser inlet curve is also built on a fixed points F and S 2 and two control points are E and the a point along S 2 M 4 line (not shown in the Figure).

Stage II
The blade angle (angle β) distribution at hub and shroud is considered for parameter study in this stage. The β angle is related to the azimuthal coordinate (angle θ) of blade curves through the following relation [1]. tan β = rdθ ds (8) in which r is the radial coordinate and s is the meridional coordinate. A five control point Bézier curve is fitted for each β angle distribution at hub and shroud ( Figure 4). To build a domain space, the slope of the control polygon at both ends are varied. The range of variation of the slopes is set to be at maximum of 5 degrees from base case. Due to β distribution variation at hub and shroud, the blade rake angle will change. The Blade Rake angle is the angle between the trailing edge and the tangential direction at impeller tip. A rake angle up to 45 • is commonly used [22] and the maximum value of rake angle is kept constant (around 54.2 • ) in this study.

Stage III
The leading edge of main blades and splitters are supposed to be changed in this stage. The main geometry of meridional profile is held constant and the leading edge starting points, H 1 and S 1 , are the lower limits for the leading edge. The upper limits for varying the position of H 1 are toward the exit by 16% of hub curve length, and for the position of S 1 are toward the exit by 10% of shroud curve length.
The leading edge position of splitter can vary ±6.5% shroud profile length and ±7.5% hub profile length related to their initial positions. Figure 5 shows the parameter bounds in this stage.

Stage IV
In this stage the Box-Behnken technique is used to create the design space based on the selected parameters in the previous stages. An optimization process can take place in this stage regarding meridional geometry, angle distribution and inlet geometry for blades simultaneously.

Geometry
A turbocharger compressor is modeled in this research. The compressor impeller has six main blades and six splitters followed by a vaneless diffuser and an overhanging volute. The original impeller and diffuser characteristics are tabulated in Table 2 and the cross section numbering is shown in Figure 6.
The splitter leading edge starts at 47% hub and 29% shroud curve length from the inlet [23]. Through evaluating the effect of each parameter on the performance, all other parameters are kept unchanged.

Modeling
The aerodynamic model has been developed using ANSYS CFX 16.0. In this model, one passage contains the inlet, impeller and diffuser and these are meshed using structured hexahedral grids. O-type and H-type schemes are used for the leading edge and trailing edge of blades respectively. The tip clearance is considered in the model and the structured grids are constructed for the gap. A boundary layer with a 50 microns thickness of the first element is implemented to maintain the wall y+ around 30. Mesh independency analysis showed that about 750,000 elements per passage is sufficient to achieve reliable results as reported. Figure 7 shows the impeller grids for the initial case at hub and shroud in which the clearance grids are detectable on the main blade tip.
At the inlet section the total pressure and temperature are set as the inlet boundary conditions to 1 bar and 300 K respectively and the mass flow rate is set at the diffuser outlet. The frozen rotor assumption is used between rotating and non-rotating mediums also explained in [24]. The shear-stress transport (SST) is used for turbulence modeling in which Reynolds stress terms are treated by k − ω model near-wall region and by k − model in the far field [25]. The energy conservation equation is solved and the solution convergence is checked by the reduction of all conservative equation residuals to the order of 1 × 10 −5 .

Effective Parameter Identification
The performance of each design point is evaluated by total pressure ratio and total to total isentropic efficiency from compressor inlet to the diffuser outlet. Pressure ratio is calculated by means of following relation: Pr = P 03 dṁ/ṁ /P 00 (9) In which the numerator is mass flow average of total pressure at the diffuser outlet. Using pressure ratio, the isentropic efficiency is calculated using Equation (10).
The normal probability plots for stage I contains 32 runs are shown in Figures 8 and 9 for the effects on η tt and Pr respectively. Among the five investigating parameters in this stage, as is expected, the impeller meridional geometry, as defined by M 1 to M 4 contains the main parameters. The diffuser inlet curvature varied by M 5 , has a negligible effect on both efficiency and pressure ratio. However the results shows that the two parameters M 1 and M 4 which adjust the hub and shroud profile at the impeller outlet respectively have the most significant effect in comparison to the others. Among the two factor interaction effects, the M 1 M 4 is apart from line in Figure 8. Both the M 1 M 4 and M 2 M 3 points in Figure 9 show that they play an important role in total pressure ratio in this case. Generally the main effects and the low order interactions are dominant and the higher order interactions are usually negligible [5].
As the M 1 to M 4 main effects are positives, this means that each parameter at high level solely improves the efficiency and pressure ratio. It should be noted that the high level for these parameters are considered in a way that the flow passage width increases. Whereas when both M 1 and M 4 at the high levels or both at the low level, efficiency and pressure ratio decrease as it has a negative value. The same analysis is true for both M 2 and M 3 on the pressure ratio.
When both M 1 and M 4 or both M 2 and M 3 are considered at high levels, it means that respectively the flow passage at the outlet or at the inlet, becomes wide and when they are at the low levels it narrows the flow passage. Therefore increasing or decreasing flow passage area by simultaneously changing the hub and shroud curves is deleterious to the performance.  In the same way, the normal probability plots for stage II, contains 16 runs, are extracted. Figures 10  and 11 are shown to investigate the effects of blade angle distribution on η tt and Pr respectively. As seen in Figure 10, B 1 , B 2 , B 3 B 4 and also all two-factor interactions with B 2 are apart from the straight line; however the other effects do not lie along the line and can be considered as significant as the mentioned ones. Whereas in Figure 10  The B 1 factor is representative for β angle distribution slope at the inlet for hub. The value of its effect on the isentropic efficiency ( Figure 10) is positive however its effect on pressure ratio ( Figure 11) is negative. This means that increasing B 1 , increases the efficiency and reduces the pressure ratio. Increasing the B 1 value is equivalent to decreasing the blade curvature as it defines the angle variation ( Figure 4). The B 2 value also shows the same effect on the pressure ratio. By considering two-factor interactions B 2 B 3 and B 3 B 4 on the efficiency, it is clear that, the efficiency will improve when the hub and shroud blade angle difference is small. In conclusion for the case II investigation, B 1 and B 2 can be considered as the most dominant factors and both efficiency and pressure ratio will improve by controlling the angle difference between hub and shroud in the defined ranges.

Effects
To better capture the effects of blade angle distribution on the fluid behavior, the streamlines for the original case and the case with {B + 1 B + 2 B + 3 B − 4 } which is shows the best efficiency are shown in Figure 12. The streamlines for one main impeller and splitter passage are shown at span 0.8 where the effect of shroud curvature which is controlled by B 3 and B 4 could be realized. The streamlines show that how the blade angle adjustment could lead to guide the flow through the passage following the blade where the fluid intensify is escaping from the tip clearance at that span. The blade outlet angle at hub and shroud for the case with B + 2 and B − 4 has the lowest difference among the test cases and it gives the most uniform flow in the span-wise direction.
Similarly, the normal probability plots for isentropic efficiency and pressure ratio are shown in Figures 13 and 14 respectively for stage III. The plots show that moving I 1 along the hub curve toward the outlet increases the pressure ratio and decreases the efficiency; however moving I 3 toward the outlet results in pressure ratio and efficiency improvement. When both I 1 and I 2 are moving simultaneously toward the inlet or outlet, it also results in the performance increasing, whereas moving I 2 towards the outlet leads to an unconventional incline at impeller inlet and reduces the pressure ratio.

Response Surface for Finding the Optimal Geometry
The above mentioned analysis for different cases show that the controlling point at the hub plays an important role in the objective functions which are the pressure ratio and isentropic efficiency. This result is may be due to following reasons:

•
The wider ranges are considered for the parameters in the hub section • The hub section variation has more effect on whole impeller geometry in comparison to the changes in the shroud section • The shroud characteristics are dependent on hub geometry hence the most effective parameters in hub also have effects on the shroud To find the optimal geometry which considers all the aforementioned analyses, just the parameters on the hub are considered for the next analysis. These are M 1 and M 2 in stage I, B 1 and B 2 in stage II and I 1 and I 3 in stage III. The other factors are set in the values which show the best performance in the previous case analysis. This leads to construction a design space for 13 (= 5 + 4 + 4) parameters. A full factorial design space with each parameter at two levels leads to 2 13 runs whereas to better capture the design space, three levels for each parameter are more helpful. However, using the full factorial technique for six parameters each at three levels consequences 3 6 runs which also imposes high computational cost.
In stage IV, as mentioned before, the Box-Benken technique is used which reduces the design space to 54 points. For each factor the levels are shown by −1, 0 and 1 in which the center value (0) is set the value was obtained by the previous results. For obtaining the desired values from the parametric study, both isentropic efficiency and total pressure ratio are considered as the objective functions by following relation [26], in which the superscript i stands for the initial value and the coefficients C 1 and C 2 are set to 0.5 for considering the same share for efficiency and pressure ratio. Minimizing f results in the best performance results in each case. Using the aforementioned objective function, in case I, M 1 M 3 M 4 run, in Case II, B 1 B 3 run and I 3 effect show the best performance results. Pressure ratio, isentropic efficiency and also the value of f for all design points are shown in Figure 16. The design points with maximum efficiency and with maximum pressure ratio are also shown in Figure 17. The results show that efficiency variation is about 1.2% and and pressure ratio variation in about 13.2% among all points. The point with maximum efficiency shows a 0.4% improvement in comparison to initial run and the point with maximum pressure ratio has 2.5% increase in this ratio.  Second-order polynomial regression is used for pressure ratio, Pr, efficiency , η, and for both simultaneously, f , using Equation (6), and the surrogate functions are noted byŷ Pr ,ŷ η andŷ f respectively. Response surface method analysis is used to find the variables {A, B, C, D, E, F} to maximizeŷ Pr andŷ η and to minimize theŷ f . Table 3 shows the results. Three extra performance analyses are performed to investigate the design parameter values {A, B, C, D, E, F} in Table 3 on the efficiency and the pressure ratio. The CFD results show an additional 0.2% improvement for pressure ratio while optimum conditions are not changed regarding the efficiency and f function. Hence related to initial condition, f 0 , the isentropic efficiency improves by about 0.4% and the pressure ratio could be increased by 2.7%. Figures 18 and 19 show the performance curve for the optimized impeller and the original impeller. Results are derived for three different rotational speeds. The design point rotational speed is shown by N and other rotational speeds are 75%N and 1.25%N. At the design point the optimized model shows 3.0% and 11.0% improvement in efficiency and pressure ratio respectively which are substantial.
The range of the optimized design is not falling off toward choke as is often the case when a new design gives greater design point performance but at the expense of range. Range improvement near surge must be determined by means of an unsteady simulation [27] of the entire compressor or by experiment; however the current results show no operating range deterioration in comparison with operating points at low and high mass flow rates.
At the lower mass flow rates the improvement in pressure ratio decreases in all rotational speeds. The isentropic efficiency curves show in lower mass flow rates, the original case has higher value in comparison to optimized case and by increasing the mass flow rate the performance of optimized impeller increases (Figure 19).   Different aspects of flow behavior can be studied to compare the original case with the optimized one. As the effect of blade angle distribution and splitter parameters on the flow streamline at different span are shown before, the pressure contour at meridional plane are investigated here. Figure 20 compares the optimized and the original case. On the left, the meridional curves are shown simultaneously as well as the leading edge positions of the main blade and splitter. The optimized geometry prepares the wider passage for flow at its mid way and it helps the fluid to diffuse more smoothly. The pressure contours at the meridional planes are shown at the design point. Beside the fact that pressure rise in optimized case is much more than the original case, it can be seen that in the original case the pressure difference between hub and shroud along the line which perpendicular to mid span is rather large. The optimized case effectively makes the pressure gradient smooth from hub to shroud as the flow passes the passage.

Conclusions
As described in this article, the geometry of meridional curves of the impeller, the angle distribution of the blade and the location of the main blade and splitter entrances were optimized. At the stage I, five parameters are considered for meridional geometry of impeller. Full Factorial analysis was performed and the most effective parameters were recognized which are the hub curve control points. The best points in this stage shows 2.5% and 7.24% improvement in efficiency and pressure ratio respectively.
The best point in the previous analysis is used was Stage II as an initial geometry to investigate the effect of the β angle distribution of the blade at hub and shroud by varying four selected parameters. Full Factorial analysis shows that the slopes at starting and at the end points of the fitted Bézier curve at hub are the effective factors. The highest pressure ratio in this stage shows 3.0% improvement and the distribution with the highest efficiency shows 0.34% increase as compared to the baseline geometry.
Considering the best performance in previous stages as an initial model in stage III, the function f (Equation (11)) is used to simultaneously consider the effect of efficiency and pressure ratio. Full factorial analysis is performed for four parameters which are the starting point locations of the main blades and splitters.
Based on previous analysis, six parameters were selected for stage IV and for decreasing the number of runs the Box-Behnken method was used to construct the design space. The efficiency and the pressure ratio for all runs were obtained and second order polynomial functions were used to find the response surfaces. Response surface analysis introduces the best value of six variables to obtain maximum efficiency, maximum pressure ratio and a minimum for f function. The results were compared to the original impeller and by considering both factors simultaneously, the optimized model shows substantial improvement in both efficiency and pressure ratio at design point condition. A very encouraging result is at different operating conditions the performance has also improved in comparison with the original case rather than be sacrificed due to improved design point performance which is often the case.