Front. Energy Power Eng. China 2008, 2(4): 541–549 DOI 10.1007/s11708-008-0077-3

RESEARCH ARTICLE

Xiaojia LIU, Fangfei NING

New response surface model and its applications in aerodynamic optimization of axial compressor blade profile

E

Higher Education Press and Springer-Verlag 2008

Abstract A parametric method for the axial compressor 2D blade profiles is proposed in which the blade geometries are defined with the parameters commonly used for blade definition, which ensures that the geometric significance is clear and an unreasonable blade profile is not generated. Several illustrations are presented to show the fitting precision of the method. A novel response surface model is proposed which regards the objective distribution function in the vicinity of a sample as normal school, and then generates the response surface function in the whole design space by a linear combination of distribution functions of all the samples. Based on this model, a numerical aerodynamic optimization platform for the axial compressor 2D blade profiles is developed, by which aerodynamic optimization of two compressor blade profiles are presented. Keywords CFD, optimization, response surface model, parametric method, compressor

1

Introduction

During the past three decades, computational fluid dynamics (CFD) played an increasingly important role in aerodynamic design of turbomachinery. As a tool that predicts flow fields and provides assessments for aerodynamic designs, CFD can help designers choose the best design from a series of candidates. But on the other hand, CFD itself cannot tell how to modify a design to achieve a better aerodynamic performance. Therefore, the design process has to be based on ad hoc practice, which is both expensive and heavily dependent on the experiences of the designers. Thus, whether the outcome achieved is at its optimum under given constraints cannot be guaranteed. Translated from Acta Aeronautica et Astronautica Sinica, 2007, 28(4): 813–820 [译自: 航空学报] Xiaojia LIU (*), Fangfei NING National Key Laboratory of Aircraft Engine, Beihang University, Beijing 100083, China E-mail: [email protected]

During the past several years, with CFD as a predictor of objective function, considerable progress has been made in optimization theories, such as shape optimization of aircrafts [1–3]. Meanwhile, much attention has been paid to the study of numerical optimization in design optimization of turbomachinery, and successful applications can be seen in literature [4–6]. Generally speaking, the goal of an optimization process is to find the minimum or maximum of an objective function within the given design space. For the optimization of a compressor blade profile, the objective function can be set as its aerodynamic property(ies), and the design space can be the geometry of the blade profile. Three key issues have direct impact on the obtained result and the efficiency of the optimization process: the method to parameterize the blade profile, the efficiency and adaptability of the optimization method, and the efficiency and accuracy of the CFD solver. This paper focuses on the two former key issues. A parametric method for the axial compressor 2D blade profiles and a novel response surface model are proposed. Combined with simulated annealing, an optimization platform was specifically constructed for the axial compressor blade profile. The optimizations for two blade profiles are presented to validate and test the efficiency of the developed platform.

2

Parametric method

A reasonable parametric method for determining blade profile should be able to represent a wide range of profiles, whilst the required parameters should be as few as possible. The representation of a wide range of profiles helps achieve the optimum whereas the fewest possible parameters is favorable for the convergence of an optimization process. Furthermore, such method should automatically be able to avoid the appearance of illogical blade profile. If this cannot be assured, (for example, if the maximal/ minimal thickness of the obtained blade profile does not match the requirement for blade strength, or if there are several inflexions of the blade profile), additional

542

constraints should be imposed which may make the objective function more complex and may even cause divergence of CFD calculations. The proposed parametric method defines a blade profile using 12 parameters: the stagger angle h, the solidity t, the maximum thickness Tm and its location along the chord Lm, the radius of the leading circle rl, the leading wedge angle dl, the leading metal angle ll, the radius of the trailing circle rt, the trailing wedge angle dt, the trailing metal angle lt, the maximum flexibility on suction surface Hs and its location along the chord Ls, as seen in Fig. 1. The suction surface is reproduced by two splines: one connects the leading edge tangential point to the maximum flexibility point on the suction surface while the other connects the maximum flexibility point to the trailing edge tangential point. At the maximum flexibility point, the continuous condition of the first-order derivative of the two splines is satisfied. At the leading/trailing edge tangential points, the slopes of the two splines are imposed by the given leading/trailing wedge angles. After the suction surface is defined, according to the maximum thickness and its location, the pressure surface can also be reproduced by the two splines, which is similar to those for suction surface, but with the middle connection point being replaced by the maximum thickness point on the pressure surface.

Xiaojia LIU, Fangfei NING

Several controlled diffusion airfoils (CDA) fitted with the proposed parametric method are shown in Fig. 2. The dashed lines denote the fitting blades, while the solid lines denote the original blades. It can be seen that this method represents the blades very well.

Fig. 2 Illustrations of several controlled diffusion airfoils fitted using parametric method (a) CDA-21 [7]; (b) CDA-MAN-GHH-1-S1 [8]; (c) CDANASA [9]

3 Optimization methodology and response surface Fig. 1 Geometrical parameters used in parametric method for axial compressor blade profile

When applying the proposed parametric method to blade profile optimization, some of the above 12 parameters can be fixed beforehand. For example, the radii of leading/trailing edge and the maximum thickness are usually determined by the requirement of blade strength, and the leading metal angle should be kept unchanged due to the direction of incoming flow and the incidence condition. It should be noted that all the parameters used in this method are those with clear geometric significance and are frequently used to describe compressor blade profiles. It is very easy to prescribe their lower and upper bounds to make sure that the represented blade profile is within a reasonable region.

An optimization strategy based on response surface model [10] is introduced. The whole optimization process is shown in Fig. 3, which includes three modules: the user interface module, the initial samples generation module, and the blade profile optimization module. In blade profile optimization module, the simulated annealing as described in Ref. [11] is used for response surface searching. A new response surface model is proposed, which is explained in detail below. Response surface is an approximation method to describe the relationship between variables or samples assigned and their corresponding objective functions. Optimization based on response surface model searches for the extremum of the response surface function rather than the realistic objective function. A response surface is usually constructed by linear or nonlinear assembly of a

Aerodynamic optimization of axial compressor blade profile

Fig. 3

543

Flow chart of numerical optimization procedure of axial compressor blade profile

series of basis functions, such as polynomial method and artificial neural networks. The former is based on the linear assembly of exponential functions, while the latter is based on the nonlinear assembly of activation functions. Normally, the larger the number of samples is, the better

the response surface converges to the real objective function (the samples should also be as scattered as possible). On the other hand, it should be realized that a group of limited samples are not able to reflect the objective function exactly. In other words, the objective function is not

544

Xiaojia LIU, Fangfei NING

clear in the design space where the samples are not covered. From this point of view, the influence area of a sample should not cover the whole design space, but just its vicinity. The commonly used polynomial method and BP networks are the methods with response surfaces constructed by global functions so that each sample could affect the response surface in the whole design space. Experience tells the authors that the optimization based on either of these two methods will easily fall into the local extremum, or even diverge. Generally, it cannot achieve the best result when there is a great difference between the basis functions to construct the response surface and the real objective function. According to the above analysis, each sample can first be assumed as one extremum of the objective function. The objective function in the vicinity of each sample can then be assumed as a distribution function. Subsequently, the response surface can be constructed by linear combination of all distribution functions corresponding to the provided samples. This new response surface method is introduced as follows. Equation (1) is used as the distribution function to fit the objective function in the vicinity of each sample: "  # L
X xl {xi:l 2 , ð1Þ fi ðxÞ~ exp { sl l~1 where x is the design variables vector; L is the number of variables; subscript ‘‘i’’ means the index of the sample, while subscript ‘‘l’’ means the index of the design variable; sl is a pre-defined constant, which defines the influence area of the sample for the given variable. The response surface can then be constructed by the linear combination of all distribution functions, as shown in Eq. (2): " 2 # N N L
X X X x {x l i:l Ci :fi ðxÞ~ Ci : exp { , ð2Þ F ðxÞ~ sl i~1 i~1 l~1 where N is the number of samples; Ci is the coefficient to be solved. With all the samples known, the linear system for Ci can be written as 3 2 2 32 3 f1 ðx1 Þ f2 ðx1 Þ . . . fN ðx1 Þ F ð x1 Þ C1 6 f1 ðx2 Þ f2 ðx2 Þ


fN ðx2 Þ 7 6 C 7 6 F ðx2 Þ 7 6 76 27 6 7 6 7:6 . 7~6 7 .. .. .. 6 .. 7 6 . 7 6 .. 7: 4 . 54 . 5 4 . 5 . . . f1 ðxN Þ f2 ðxN Þ




fN ðxN Þ

CN

F ðxN Þ

After Ci is solved, the response surface function can be obtained as defined by Eq. (2). In this model, the predefined sl could have a different value for each variable. Assuming that the range of a specific design variable is 0 to 1, generally, sl should have a slightly larger value if the objective function has just one extremum about this design variable, for example, sl 5 1–2;

otherwise, it should be assigned as a smaller value, for example, sl 5 0.2–0.5. In actual applications, the value of sl should be decided by numerical experiment. This response surface model is advantageous in that it can reproduce arbitrary objective function; and with a new sample produced, the response surface can be improved in the vicinity of this sample. On one hand, the convergence of the optimization can be accelerated; on the other hand, the design space without samples can be explored. An illustration shows the efficiency of this response surface model. The model is used in the optimization process to search the universal minimum value of the function with a single variable as follows:   y~1z expð{1:9xÞ: sin 10x2 , 0fxf1: First, 3 samples are generated randomly. Based on these samples, the response surface is constructed. Second, the minimum point of the response surface is found using simulated annealing. The real objective function of this point is calculated, and a new sample is obtained which is then added to the sample set to construct a new response surface for the next round of search. The second step is repeated until the minimum point does not change any more. At this time, the minimum point of response surface is equal to the minimum objective in reality. In this case, the optimization converges after 8 circles, as seen in Fig. 4. The convergence history shows that the number of initial samples can be very small. The characteristic of this model makes the minimum point tend to appear at the region far from the existing samples, which keeps the distribution of the samples much more uniform and benefits global searching.

4 Numerical simulation A self-developed Multi-block Aerodynamic Prediction (MAP) code based on finite volume method for spatial discretization is used as the flow solver. The convective flux is evaluated with AUSMD/V method, while the viscous flux is calculated with traditional centre differentiation. GMRES algorithm is used to solve the implicit time discretized governing equations. Moreover, AGS transient model is adopted to simulate the turbulent transient flow. More details about numerical methods can be referred to in Refs. [12–14]. In Fig. 5, blade surface isentropic Mach number distributions of a 2D compressor cascade at different Reynolds numbers, showing the accuracy of the CFD code. The chord length of this cascade is 100 mm; the stagger angle is 28.4u; the solidity is 1.63. H-O-H grid (the O-type around the blade and Htype at the up/down stream area) is used and the total number of the grid is 5500. The turbulent intensity of the inlet flow is 5.6%; the incidence corrected with the

Aerodynamic optimization of axial compressor blade profile

Fig. 4

545

Process of searching for global optimum using proposed response surface model

experiment is 1.8u; the inlet Mach numbers are 0.700 and 0.775, respectively; the chord Reynolds numbers are 5.5 6 105 and 1.9 6 105, respectively. The boundary layer remains attached along the suction surface at high Reynolds number, while the separation bubble appears

on the suction surface at low Reynolds number. In both cases, the CFD results match with the experiments very well, especially at the low Reynolds number where the CFD accurately predicts the presence of the separation bubble, as seen in Fig. 5(b).

546

Xiaojia LIU, Fangfei NING

Fig. 5

Blade surface isentropic Mach number distributions of compressor cascade at different Reynolds numbers (a) Re 5 5.5 6 105; (b) Re 5 1.9 6 105

5 Optimization of axial compressor blade profiles The optimization of two axial compressor blade profiles with the developed optimization platform is presented in this section. One blade profile is the midspan of a stator and the other is the midspan of a rotor. Both profiles are obtained after deliberate design. The work presented here should further improve their performances. The objective function is defined as follows: OBJ 5 optimization objective function + penalty function, where, the optimization objective function is set as the total pressure loss coefficient, while the penalty function is the difference between the absolute outlet flow angles of the new blade profile and the initial blade profile. Therefore, the total optimization objective should achieve minimum loss while keeping the outlet flow angle unchanged. The parametric representation method used in the optimizations is that introduced in Section 2, and the design variables and their corresponding ranges can be found in Table 1. The other geometric parameters, such as the stagger angle, the leading edge and trailing edge radii, the maximum thickness, the leading edge metal angle, and the solidity, remain unchanged. The stagger angle is decided by the multistage matching condition and the Table 1 List of variable design parameters of blade profile and their variation extensions design variable

range(relative to initial blade)

trailing edge metal angle/(u) leading edge wedge angle/(u) trailing edge wedge angle/(u) suction surface maximum flexibility

¡10 ¡5 ¡5 ¡20% initial maximum flexibility ¡10% chord

relative location of suction surface maximum flexibility relative location of maximum thickness

¡10% chord

blade stacking. The leading/trailing edge radii and the maximum thickness are fixed according to the strength requirement. The leading metal angle is designed according to the inlet flow angle and the incidence angle under the 3D multi-row environment. The solidity which relates to the chord length of the blade profile and the number of the blades is not considered here. Before the optimization, a set of initial samples has to be generated. In the present work, 20 initial samples are selected using the method introduced as follows to ensure that the sample distribution is as uniform as possible in the search space. Assuming that k samples have been generated, then the remaining samples can be generated as follows: (1) select design variables randomly. (2) if (min|xk+1 – xi| . d (i 5 1,2,???,k)), then the design variables is accepted as a new sample; else go to (1). end where, x is the design variable; | ? | is the Euclid norm; d is the constant assigned beforehand, by which the distances among the samples can be controlled. The value of d can be assigned by numerical experiments. 5.1

Optimization result of stator

The main geometric parameters of the initial stator are: a chord of 24 mm; a blade angle of 46.7u; a maximum thickness that is 7.2% of the chord, and a location that is 49% of the chord; a stagger angle of 28.6u; a local solidity of 1.32. The main aerodynamic parameters are: the inlet Mach number 0.7; the design incidence angle 0u; the diffusion factor 0.53. The convergence history of the optimization is shown in Fig. 6. The two lines in Fig. 6 represent the optimum objective function provided by the response surface and the actual value obtained at the optimum point by the CFD, respectively. The final optimum can be obtained when the two calculations become equal to each other. In this case, the optimization converges after 15

Aerodynamic optimization of axial compressor blade profile

547

distributions of boundary layers along the suction surfaces of the original and the optimized stator blade profiles are shown in Fig. 8. Though the separation is not totally removed, the separation point of the optimized blade profile moves downstream. The width of the wake downstream of the optimized blade profile is much smaller than that of the initial one, as seen in Fig. 9. The outlet flow angles of the initial profile and the optimized profile remain the same, which are 18.46u and 18.49u respectively. The reduction of the loss coefficient is about 20%, which is from 0.032 to 0.026 with the given constraint conditions. 5.2 Fig. 6

Optimization history of stator blade profile

iterations. The Mach number distributions of both the initial blade profile and the optimized blade profile are shown in Fig. 7. It can be seen that the suction surface maximum flexibility of the optimized profile is reduced and the blade loading moves upstream. The shape factor

Fig. 7 Comparison of original and optimized stator blade profiles and blade surface isentropic Mach number distributions

Fig. 9

Optimization result of rotor

The main geometric parameters of the initial rotor are: a chord of 38 mm; a blade angle of 38.7u; a maximum thickness that is 7.0% of the chord, and a location at 50% of the chord; a stagger angle of 40.2u; a local solidity of 1.53. The main aerodynamic parameters are: an inlet relative Mach number of 0.7; a design incidence angle of 21.0u; a

Fig. 8 Shape factor distributions of boundary layers along suction surfaces of original and optimized stator blade profiles

Mach number distributions of original and optimized stator flow fields (a) Original design; (b) optimized design

548

Xiaojia LIU, Fangfei NING

diffusion factor of 0.56. The optimization converges after 40 iterations, as seen in Fig. 10. The blade surface isentropic Mach number distribution in Fig. 11 shows that there is a local supersonic area on suction surface and the blade loading moves upstream. The momentum thickness distribution of boundary layers along the suction surfaces of the original and the optimized rotor blade profiles are shown in Fig. 12. It can be seen that the boundary layer thickness leaving the solid surface is considerably reduced. The Mach number distributions of the initial blade profile and the optimized profile are shown in Fig. 13. Compared with that of the initial profile, the separation near the trailing edge of suction surface and the width of the downstream wake are obviously subdued in the optimized profile. The outlet flow angles, which are 50.49u for the initial blade profile and 50.29u for the optimized profile, respectively, remain the same, while the loss coefficient is reduced by about 16%, from 0.086 to 0.072.

Fig. 10

History of optimization process of rotor blade profile

Fig. 13

Fig. 11 Comparison of original and optimized rotor blade profiles and blade surface isentropic Mach number distributions

Fig. 12 Momentum thickness distributions of boundary layers along suction surfaces of original and optimized stator blade profiles

Mach number distributions of original and optimized rotor flow fields (a) Original design; (b) optimized design

Aerodynamic optimization of axial compressor blade profile

6

Conclusions

1) A parametric method for the axial compressor 2D blade profile is proposed. Several illustrations show that it could describe different blades with a few free parameters, and that it is particularly able to avoid the appearance of an unreasonable profile. If a parametric method is expected to automatically avoid the production of unreasonable blade profile, the parameters with obvious geometrical significance for blade profiles should be used, and a low order curve should be adopted to avoid unexpected inflexion. 2) A novel response surface model is proposed in which the attribute of each sample is limited to its vicinity. The characteristics of this model are as follows: a. It doesnot require many initial samples. The new sample tends to appear in the search space with few samples so that the response surface can be updated efficiently. b. Theoretically, it is able to fit in the arbitrary objective function. c. According to the authors’ experience, it is insensitive to error when evaluating the actual objective function, i.e., the calculation error of actual objective function would not affect the response surface in the whole search space. It just affects that within its own vicinity. Numerical optimization of the two axial compressor blade profiles proves that this response surface optimization methodology is efficient. It can be seen as a supplement to the existing optimization methods.

References 1. Jameson A, Martinelli L, Alonso J J, et al. Simulation based aerodynamic design. IEEE Aerospace Conference, Big Sky, MO, 2000

549 2. Zhu Ziqiang, Fu Hongyan, Yu Xinri, et al. Multi-objective optimization design of airfoil and wing. Science in China, Series E, 2003, 33(11): 999–1006 (in Chinese) 3. Xiong Juntao, Qiao Zhide, Han Zhonghua. Optimum aerodynamic design of transonic wing based on response surface methodology. Acta Aeronautica Et Astronautica Sinica, 2006, 27(3): 399–402 (in Chinese) 4. Pierret S, Van den Braembussche R A. Turbomachinery blade design using a Navier-Stokes solver and artificial neural network. ASME Paper 98-GT-4, 1998 5. Eyi S, Lee K D. Turbomachinery blade design via optimization. AIAA Paper 2000–0740, 2000 6. Koller U, Monig R, Kusters B, et al. Development of advanced compressor airfoils for heavy-duty gas turbines, part I: design and optimization. ASME Paper 99-GT-95, 1999 7. Saha U K, Roy B. Experimental investigations on tandem compressor cascade performance at low speeds. Experimental Thermal and Fluid Science, 1997, 14(3): 263–276 8. Steinert W, Eisenberg B, Starken H. Design and testing of a controlled diffusion airfoil cascade for industrial axial flow compressor application. ASME 90-GT-140, 1990 9. Ho Y K, Walker G J, Stow P. Boundary layer and NavierStokes analysis of a NASA controlled-diffusion compressor blade. ASME 90-GT-236, 1990 10. Shyy W, Papila N, Vaidyanathan R, et al. Global design optimization for aerodynamics and rocket propulsion components. Progress in Aerospace Sciences, 2001, 37: 59–118 11. Tsallis C, Stariolo D A. Generalized simulated annealing. Physica A, 1996, 233: 395–406 12. Ning Fangfei. Numerical investigations of flows in transonic compressors with real geometrical complexities. Dissertation for the Doctoral Degree. Beijing: Beijing University of Aeronautics & Astronautics, 2002 (in Chinese) 13. Ning F F, Xu L P. Numerical investigation of transonic compressor rotor flow using an implicit 3D flow solver with oneequation Spalart-Allmaras turbulence model. ASME Paper 2001-GT-0359, 2001 14. Abu-Ghannam B J, Shaw R. Natural transition of boundary layers - the effect of turbulence, pressure gradient, and flow history. Journal of Mechanical Engineering and Science, 1980, 22(5): 213–228