Jinxin CHENG,Jiang CHEN,Hang XIANG

School of Energy and Power Engineering,Beihang University,Beijing 100083,China

KEYWORDS Aerodynamic optimization;Bezier surface;Compressor;Global optimization;Surface parametric control

Abstract An aerodynamic optimization method for axial flow compressor blades available for engineering is developed in this paper.Bezier surface is adopted as parameterization method to control the suction surface of the blades,which brings the following advantages:(A)significantly reducing design variables;(B)easy to ensure the mechanical strength of rotating blades;(C)better physical understanding;(D)easy to achieve smooth surface.The Improved Artificial Bee Colony(IABC)algorithm,which significantly increases the convergence speed and global optimization ability,is adopted to find the optimal result.A new engineering optimization tool is constructed by combining the surface parametric control method,the IABC algorithm,with a verified Computational Fluid Dynamics(CFD)simulation method,and it has been successfully applied in the aerodynamic optimization for a single-row transonic rotor(Rotor 37)and a single-stage transonic axial flow compressor(Stage 35).With the constraint that the relative change in the flow rate is less than 0.5%and the total pressure ratio does not decrease,within the acceptable time in engineering,the adiabatic efficiency of Rotor 37 at design point increases by 1.02%,while its surge margin 0.84%,and the adiabatic efficiency of Stage 35 0.54%,while its surge margin 1.11%after optimization,to verify the effectiveness and potential in engineering of this new tool for optimization of axial compressor blade.

1.Introduction

The aerodynamic design of axial flow compressors is the first key step in the design of aeroengines.The goal of aerodynamic design is to achieve‘‘three high performance”of‘‘high efficiency,a high pressure ratio,and a high surge margin”,but the‘‘three high performance”is often interrelated and contradictory.The traditional aerodynamic design method of axial flow compressors is to combine one-dimensional design and analysis,two-dimensional design and analysis,quasi threedimensional design,and three-dimensional analysis,and its disadvantage is that it highly relies on expert experience and can only be compared and improved in a limited number of scenarios.In order to improve compressor design performance with as little manual intervention as possible,optimization methods have been introduced into aerodynamic design of compressors since the 1980s.

The key to the aerodynamic optimization of axial flow compressor is twofold:one is geometrical parameterization;the other is the optimization algorithm.Geometrical parameterization is a method for controlling the geometrical deformation of a blade or a flow channel according to several parameters in an optimization process.The traditional method of compressor blade parameterization is based on changes to radial sections. One change involves the geometry of each section being held constant while the relative position changes;that is,the axial and circumferential offsets of the section constitute the bow and sweep of the blade.1-3The other approach is to change the shape of central arced curve and thickness distribution in each section,4-7or to change the suction or pressure side profile of each section by means of free curves,8-10so as to change the blade geometry. The disadvantages of the method mentioned above include a large number of optimization parameters and non-smooth surface.In 2003,Burguburu and le Pape11used a Bezier surface to parameterize a singlerow transonic axial flow rotor.This surface parametric method has the advantages including fewer design parameters,better physical understanding,smooth surface and easy to guarantee the mechanical strength,and it is an important development direction in blade parameterization currently.In the aerodynamic design of axial flow compressors,aerodynamic optimization algorithms have developed from local optimization to global one.In 1983,Sanger12combined optimization technology with compressor aerodynamic designs for the first time,using the general purpose control algorithm(a local optimization algorithm),which began the era of aerodynamic compressor optimization design.In 1999,Koller and Monog13used an optimization algorithm,which combined local optimization with global one,to optimize the subsonic profile,increasing the incidence range of the compressor blades.After 2000,Ashihara14and Oyama15et al.adopted genetic algorithms to optimize three-dimensional(3D)compressor blades,improving the compressor's efficiency and surge margin at its design point.During the last two decades,artificial intelligence algorithms,as a global optimization algorithm,have been widely used in aerospace and other engineering fields.New artificial intelligence algorithms, including the artificial ant colony algorithm16(proposed in 1992), particle swarm algorithm17(proposed in 1995)and Artificial Bee Colony(ABC)algorithm18(proposed in 2005), have developed rapidly. The ABC algorithm has a unique advantage in terms of searching ability among the intelligent algorithm.19It performs both global and local searches in each iteration,combining the advantages of global optimization(accuracy)and rapid convergence(speed),making it suitable for constrained,multi-parameter,and multi-objective global optimization of complex systems,such as the aerodynamic design of multi-stage axial flow compressors.However,the ABC algorithm does not make full use of the overall food source information in its search process,and it can be further improved in terms of search capability and convergence speed.

3D aerodynamic optimization of axial flow compressor is a typically High-dimensional, time-consuming-Expensive and Black-box(HEB)problem.The core of the engineering optimization for HEB problems is to adopt kinds of strategies to achieve the relatively better optimal result,20at expense of the precision of optimal solution,within limited time by using the existing computing resources.Currently used optimization strategies include decomposition method,21-23surrogate model method,24-26non-gradient optimization method27-28and direct optimization method.29-31Ref.23points out that no optimization strategy is completely superior,and different combination optimization strategies should be adopted for different situations.

This study develops a new tool for engineering optimization of axial compressor blades. The tool adopts the efficient Improved ABC(IABC)algorithm to operate finite number of iterations of bee colony (non-gradient optimization method),and utilize the Bezier surface parametric control method to modify the suction surface of blades(decomposition method),which significantly reduced design variables,and employs the strategy of using coarse mesh in the optimization process whilst fine mesh to calculate the optimal solution(direct optimization strategy), combining with the verified CFD simulation method,to form a new engineering global optimization method.This method is applied to optimizing the aerodynamic performance of a single-row transonic rotor(Rotor 37)and a single-stage transonic axial flow compressor(Stage 35)in order to verify the effectiveness,exploring the acceptable engineering optimization approach for multi-stage axial flow compressor with constraint and multi-parameter.

2.Surface parameterization control method

The essence of the surface parametric control method is to add a Bezier surface on the original one,using the control points of the Bezier surface as the optimization variables to control the blade geometry in the optimization cycle,as shown in Fig.1.In Fig.1,LE and TE refer to the leading edge and trailing edge.During the adding process,the four vertices of the Bezier surface correspond to the four vertices of the original one.The high-order continuity of each point on the Bezier surface32ensures the smoothness of the optimized surface,which can reduce the flow loss caused by the partial roughness of the blade.

Fig.1 Surface parameterization diagram.

Fig.2 Bezier surface parameterization process.

As shown in Fig.2,the Bezier surface parameterization process can be divided into the following three steps:

Step 1.To parameterize arc length for each section of the original surface.Since the Bezier surface is a unit surface in the computational domain, the original blade profile that needs to be optimized is first transformed into the unit plane in the calculation domain through arc length parameterization,so that points in the physical and computational domains can be mapped one by one. Arc length parameterization is achieved by

Step 2.The Bezier surface function is used to calculate the variation of each point in the computation domain of the original surface.The Bezier surface function is defined as

In Eq.(3),R represents the distance of movement of each point in the computational plane that is perpendicular to the direction of the computational plane;Pk,lrefers to the control point of the Bezier surface,and the number of control points is(m +1)× (n +1);are Bernstein function determined by Eq.(4),where v and u are the two coordinate axes of the Bezier surface in computational domain,respectively,and their range of variation is[0,1];is the number of combinations determined by Eq.(5).

Step 3.The R value is combined with the direction of movement of each point of the original blade profile to obtain a new one.

3.IABC algorithm

3.1.Basic principles of ABC algorithms

The ABC algorithm is a global optimization algorithm,which is inspired by the process of bees looking for food source,as shown in Fig.3.

The bee colony involved in the optimization process consists of three kinds:employed bees,onlooker bees and detector bees.In the initial phase,half of the bee population are employed bees and the other half are onlooker bees.The locations of food sources collected by employed bees correspond to a set of feasible solutions,and the food source concentration represents the fitness of a feasible solution.The process of the ABC algorithm is as follows:

Step 1.Initialize the food source.First,a feasible solution of N D-dimensional vectors Xi(i=1,2,···,N)is randomly generated,where D is the number of optimization variables.The formula for randomly generating the initial food source is as follows:

where j ∈ {1 ,2 ,···,D}stands for a component of the Ddimensional solution vector;represents the jth component of the ith individual;Rand(0,1)is a random number between(0,1);anddenote the smallest and largest feasible solutions of the jth component of the N feasible solutions,respectively.

Step 2.Calculate the initial food source concentration;that is,the fitness of the feasible solution of the group,record the best value,and sort according to the fitness value.The bees whose fitness value of the food sources rank the first half serve as employed bees,while the remaining bees serve as onlookers.

Step 3.Each bee adopts Eq.(7)to explore other food sources near the original food source,and update the fitness value.If the new food source has higher fitness value,it will replace the original one.If it is lower,the bee will explore new food source next time.

Step 4.In accordance with the Russian Roulette Law,the onlooker bees select a food source with probability proportional to the fitness of food source explored by employed bees,then the onlooker bees explore new food sources in the vicinity of the existing sources.If the new food source has a higher fitness value,onlooker bees become employed bees and replace the original food source with the new one.If the opposite is true,onlooker bees do the next exploration.The onlooker bees select the food source with the probability calculated by

Fig.3 Flowchart of a standard ABC algorithm.

where Nestands for the number of employed bees,that is,is the probability of choosing a food source;f(Xi)and f(Xm)stand for the fitness of the ith and mth individuals,respectively.

Step 5.If the number of explorations by employed bees or onlooker bees exceeds a certain limit and no food source with higher fitness is found,the employed bee or the onlooker bee will be transferred to a detector bee,and the detector bee will regenerate a new food source according to Eq.(6).After the new food source is generated,the detector bees are converted to employed bees.

Step 6.Record the best food source so far,and skip to Step 2 until the out-of-loop condition is met,and output the best food source location.

3.2.IABC algorithm principles

In the ABC algorithm,both onlooker and employed bees use Eq.(7)to search for new food sources.However,this formula only uses local random information of the total bee colony,which leads to the lack of capabilities in global exploration and exploitation.Eq.(7)is modified in the IABC algorithm,employed bees adopting Eq.(9)which contains the information of global optimum point to search for new food sources whilst onlooker bees adopting Eq.(10)which contains the information of local optimum point to search for new food source.Better utilizing the comprehensive information of the last generation by Eqs.(9)and(10)is the core cause for improving performance of the IABC algorithm.

where j ∈ {1 ,2 ,···,D},and d(i ,t)represents the Chebyshev distances between Xiand Xt.

The neighborhood of Xiis determined by Eq.(12),and the formula shows that when the Chebyshev distance between Xiand Xtis less than the product of the neighborhood radius r and the average Chebyshev distance,Xtis in the Xineighborhood,otherwise it is not.

where mdiis the average Chebyshev distance between Xiand the total onlooker bee colony,S is the area adjacent to Xi,r is the radius of the adjacent area,and experience shows that when r=1,the algorithm has the best convergence.

where fit(·)is the individual fitness of each bee.

3.3.IABC algorithm verification

To verify the superiority of the IABC algorithm compared to general Genetic Algorithm(GA)and ABC algorithm,two benchmark functions(sphere and Griewank functions)are utilized.Both benchmark functions have a variable dimension of 30,a function minimum of 0,a sphere function range of[-5.12,5.12],and a Griewank function variable of[-600.0,600.0].We set the number of bee colonies for the two functions to 200,the number of iterations to 200,and the maximum number of exploration to 100.Table 1 and Fig.4 show the result that both the convergence speed and global optimization ability of IABC algorithm have a significant increase compared to GA and ABC algorithm.

4.Engineering global optimization for Rotor 37 and Stage 35 based on an IABC algorithm and surface parametric control method

In order to verify effectiveness and superiority of the surface parametric control method and the IABC optimization algorithm in the aerodynamic optimization of an axial flow compressor,a global optimization platform combining surface parametric control method,the IABC optimization algorithm with 3D CFD numerical simulation was constructed,firstly applying it on Rotor 37,a single-row transonic rotor,and secondly Stage 35,a single-stage transonic compressor as the optimization object for verifying.

Rotor 37,with design speed of 17188.7 r/min and 36 blades,is one of four high-pressure compressor inlet stages designed by NASA's Glenn Research Center in the 1970s,and it has geometric data and detailed experimental data that can be found in Ref.33

Stage 35,with design speed of 17188.7 r/min,36 rotor blades and 46 stator blades,is a low-aspect-ratio single-stage transonic axial flow compressor developed by NASA's Glenn Research Center in 1978,and the compressor has detailed geometric and experimental data available in Ref.33

4.1.Optimization method

4.1.1.Blade parameterization

Optimum blade obtained in traditional way controlling the stacking line has characteristic of bow and sweep.A large torque will be produced for rotors when rotating at high speed,difficult to guarantee the mechanical strength of the blade.34Bezier surface parametric control method hardly change the characteristic of bow and sweep for blades,leading to overcoming this shortcoming.

Table 1 Comparison of test results among three algorithms.

Fig.4 Comparison of performance of ABC and IABC algorithms.

The key problem of aerodynamic optimization of compressors is to solve the contradiction between design space and time-consumption.The design space increases exponentially with the increase of the design variables,even suffering from‘‘dimension curse”,35leading to the failure of optimization.So the requirement in the practical engineering is to minimize the number of design variables under the condition of guaranteeing the sufficient design space.

Under normal circumstances,compared to the pressure surface of the blade,the loss caused by airflow on the suction surface is the main source of efficiency reduction.Therefore,only the suction surface of the blade is selected as the optimized surface in the optimization cycle to achieve the purpose of dimensionality reduction.The blade suction surface is parameterized using a 6×3 order Bezier surface.Considering that the geometry near leading and trailing edge has a non-negligible influence on flow field,as shown in Fig.5,seven control points are set in the ξ direction,with positions at the leading and trailing edges(0%and 100%)and at the 10%,30%,50%,70%and 90%positions whist four points are set in the η direction,with position at the 0%,20%,50%and 100%.Due to the guarantee that the first derivative of the connection point,which is between the suction side and pressure side at the leading and trailing edge,is continuous,the first two points(ξ1,ξ2)and the last two points(ξ6,ξ7)of each radial height are set as the fixed points.Taking into account the mechanical strength of the rotor components,the points at the 0%radial height is set as fixed points.The black and red points indicate the fixed point and the active point in the optimization process,respectively,so only nine points per blade are selected for optimization,reducing the computing costs significantly.

4.1.2.Optimization algorithm

Fig.5 Distribution of control points on Bezier surface.

To reduce computing costs of 3D aerodynamic optimization for compressors in terms of algorithm,the traditional method is to adopt adjoint algorithm and surrogate model.Adjoint algorithm was introduced to aerospace field by Jemeson in 198836,37and then this method was also applied in the optimization of turbomachinery,obtaining lots of scientific research result.38-40The advantage of adjoint method is that the time spent on optimization is independent of the number of design variables,but the disadvantage is only local optimization.The essence of surrogate model is to obtain approximate function of the time-consuming flow field simulation by sample training,reducing the computing cost significantly.Two disadvantages of the surrogate model are as follows:one is that when the design variables gradually increase,the number of samples trained to obtain approximate functions increases rapidly,which increases the calculation cost and even exceeds the scope allowed by the engineering;the other one is that in the refined optimization of the original blade with small performance improvement space,the error result is likely to occur due to the limited accuracy of the approximation function.

The new optimization method proposed in this paper adopts IABC algorithm for optimization.The number of bees in the algorithm can increase correspondingly with the increase of optimization variables.The number of iteration steps can be set artificially and depends on the actual needs.In this way,the optimized solution(not the optimal solution)can be found globally within the engineering allowed time cost,which overcomes the disadvantages of adjoint algorithm and surrogate model algorithm.Considering the time cost,based on the existing experience,this optimization sets the colony to iterate for 3 times to obtain the relative optimized solution.

4.1.3.Optimization objective function and constraints

In order to save computing cost,according to experience,we set the point which has relative high back pressure as the optimization condition to achieve the goals that both the adiabatic efficiency and the surge margin could improve.The objective function in the practical optimization process was set as follows:

max f=eff

The constraint is as follows:

where eff means the adiabatic efficiency;TPR and TPRoriare the total pressure ratio and the original total pressure ratio in the optimization process, respectively; mass and massoriare the flow rate and the original flow rate in the optimization process,respectively;minus is a very small value which is set artificially;xiis the optimization variables;andare the upper and lower limit of optimization variables,respectively.

We set the adiabatic efficiency under the optimization condition as the optimization goal,and set a relative pressure ratio change of not more than 0.5%and the flow rate not decrease as the strong constraint condition.In the optimization process,we set the objective function to the minus value when the solution does not match the constraints,thus eliminating this position of food source.For the purpose of increasing the number of feasible solutions, the exploring times of bees can be increased properly when the strong constraint occurs.

4.1.4.Optimization process

Fig.6 Flowchart of optimization process.

The optimization process is shown in Fig.6.Firstly,the program reads the geometric optimization variables of the blade,then initializes the food source,obtains the initial solution,and then calculates the fitness of each initial solution.The calculation of fitness consists of blade geometry generation,grid generation and flow field calculation.At this point,if the fitness reaches the condition of an exit loop,the optimization is finished and the optimized blade geometry is output.Otherwise,the IABC algorithm is used for optimization exploration,so as to give a new feasible solution and complete an iteration.The loop continues until the exit condition is met(i.e.,it converges or reaches the maximum number of iterations),and thus the final optimal blade is obtained.

4.2.Numerical simulation

4.2.1.Numerical methods

Flow field calculation is the premise of performance evaluation in the optimization process.The numerical calculation tool is the Fine Turbo module in NUMECA_V9.1,and Spalart-Allmaras(S-A)one-equation model is selected as the turbulent model.The fourth-order explicit Runge-Kutta method is used to solve the equation for time advancement,and the central difference scheme and artificial viscosity are used to control the false oscillations near the shockwaves and eliminate other minor oscillations.The convergence is accelerated by a local time-step multi-grid technique and implicit residual error method.

The inlet boundary conditions are: total pressure 101325 Pa,total temperature 293.15 K,and incoming flow axial direction.The solid wall is a no-slip boundary condition.The outlet boundary conditions are the given back-pressure,and the back-pressure at the outlet is gradually adjusted from the jam point to the near stall point during the calculation process.The near stall point is the highest back-pressure point before the divergence,and the blockage point is the minimum back-pressure point before the divergence is calculated.

4.2.2.Numerical method check

Stage 35 with detailed aerodynamic experimental data is adopted at design speed to verify numerical method.As shown in Fig.7(a),the experimental highest pressure ratio is 3.7%larger than the simulated one,whilst the blockage mass flow is 1.4%smaller.Fig.7(b)shows that the experimental highest efficiency is 1.1%greater than the simulated one,whilst the surge margin is 8.44%higher.According to the analysis,although there is a relative error between the CFD calculation and the experiment, the curves obtained by these two approaches have the same trend,and relatively good calculation accuracy is obtained,which ensures the reliability of the flow field calculation in the aerodynamic optimization cycle.

4.2.3.Grid generation and grid independence verification

Stage 35 is also used for grid independence verification.The grids are generated automatically by the Autogrid5 module in NUMECA_V9.1.We set the first layer near-wall grid spacing as 0.001 mm,to guarantee the near-wall grid y+≤5.The tip clearance of the rotor blade was set as 0.4 mm,the hub clearance of the stator blade as 0.4 mm,and the grid topology as 4HO.

Four sets of grids are adopted for the Stage 35 reference blade,and all grid qualities meet the calculation requirements.A comparison of the flow field calculation results of the four sets of grids is shown in Fig.8.The grid numbers of the rotor blade with the four sets of grid templates(Mesh 1,Mesh 2,Mesh 3 and Mesh 4) are 320000, 680000, 1020000 and 1840000,and the stator blade 340000,750000,1100000 and 1880000.The difference between the flow field calculation results of the third and fourth sets of grids is very small,so the third set of grids can meet the requirement of grid independence.

For the three-dimensional optimization of the compressor,in order to ensure a certain calculation accuracy in the optimization process,and at the same time save the computational cost,this paper uses the‘‘rough grid optimization,fine grid verification”method used in Ref.41:first using the second set of grids to optimize,and then using the third set of grids for flow field analysis.This method saves about 1/3 of the time.

4.3.Rotor 37 optimization verification

Rotor 37 is optimized by the new optimization method,and nine design variables are set.According to experience,the design space of each variable is specified to be[-0.05,0.1],where the negative values represent the direction of expansion at the suction surface of the blade,while the positive values are the direction of adduction at the suction surface of the blade.The Rotor 37 blade tip clearance is set to 0.36 mm and the grid number to 680000.The flow field will be calculated after obtaining the optimized blade using 1.02 million grids.We set the scale of the bee colony in IABC algorithm to 50,the maximum number of iteration to 3,and the maximum number of exploitation to 3.

Fig.7 Comparison of performance of Stage 35 between experiment and simulation.

Fig.8 Grid independence verification.

The computer used to run the optimization process has an Intel Core i7 3.07 GHz processor and 2 GB RAM.In the case of parallel computation using 5 CPUs,each generation of the bee colony optimization takes 8 h and a total of three optimization iterations are performed.Hence,the total runtime is 24 h,within the acceptable range of engineering.The history of optimization cycle is shown in Fig.9,and fopt-forimeans the difference between the optimal efficiency and the original efficiency after each iteration.

4.3.1.Comparison and analysis of optimization results

As shown in Table 2,at the design point and design speed,after optimization,the flow is increased by 1.82%,the adiabatic efficiency by 1.02%,and the surge margin by 0.84%whilst keeping the total pressure ratio almost constant.

Fig.10 shows a comparison of the calculated aerodynamic performance of Rotor 37 before and after optimization at the design speed.From Fig.10(a),it can be observed that under the same pressure ratio,the flow of the optimized blade is larger than that of the original blade in the full flow range,with the maximum difference not exceeding 2.0%,within the acceptable range of engineering.In Fig.10(b),the adiabatic efficiency of the optimized blade is higher than that of the original blade in the full flow range,and the maximum difference is 1.2%.

Fig.9 History of optimization cycle of Rotor 37.

Fig.10 Comparison of aerodynamic performance before and after optimization of Rotor 37.

As shown in Fig.11,compared with the original blade,the suction surface of the optimized blade has an outward expansion at the trailing edge,and with the increase of the radial height,the extent of the outward expansion area in the chord direction increases,and the thickness of the blade in outer expanding area increases,while in other areas,the suction surface of optimized blade is adducted and the blade thickness is reduced.

Fig.12 shows the variation of each point of the optimized blade relative to the original one on the computational plane,where the positive value represents the adduction direction and the negative value represents the outward expansion direction.The intersection area of the range above 50%radial height and 20%-50%range of chord length direction has the largest adduction amplitude of the blade.Since the geometric variation of the blade's suction surface is very small,the amount of change before and after optimization is multiplied by 10,as shown in Fig.13,so that the difference of geometry of the section at different heights(h/H)can be seen more clearly.The optimized airfoil is adducted at the hub section and is‘‘S-shaped”at the middle and tip sections,whose inflection point at the middle and tip section is located at about 70%of the chord length from the leading edge,and the other area is adducted except for the expansion near the trailing edge.

Fig.11 Comparison of geometry of suction surface before and after optimization of Rotor 37.

4.3.2.Comparison and analysis of flow field before and after optimization

Fig.12 Distribution of changes on suction surface of Rotor 37.

Fig.13 Comparison of geometry of hub,middle and tip sections before and after optimization of Rotor 37(multiply magnitude of change by 10).

Fig.14 Comparison of incidence angle distribution of incoming airflow in radial direction before and after optimization of Rotor 37.

The flow field before and after optimization is analyzed in combination with Figs.14-16.As the positive incidence angles of the incoming airflow at each section in the radial direction of the optimized blade are all smaller than the original one(Fig.14),the accelerated distance of the airflow at the leading edge of the optimized blade shortens,and the intensity of oblique shockwave near the leading edge decreases,and thus the loss of shockwave and that caused by interaction between the shockwave and the boundary layer decrease.According to the analysis in Figs.13,15(a),(b)and 16(a),the suction surface of the optimized blade at the hub section is slightly adducted,which slows down the airflow acceleration at the suction side,slightly decreasing the intensity of the channel shockwave,thus decreasing the loss of shockwave slightly.And the deceleration of airflow acceleration also pushes back the channel shockwave position slightly,reducing the range of airflow separation caused by the boundary layer after the shockwave,and thus the separation loss near the trailing edge of the hub area of the optimized blade is slightly reduced.As shown in Figs.13,15(c)-(f)and 16(b)-(c),compared to the hub section, the middle and tip one have more obvious adducted profiles at the suction side,which locate ahead of 70%chord length from leading edge,and thus the reduction of the loss,including the shockwave loss caused by the deceleration of airflow acceleration,the loss caused by the interaction between the shockwave and boundary layer and the separation loss,are more obvious.

Due to the reduction of separation areas of the boundary layer and the degree of airflow separation,which are at the hub,middle,and tip sections of the optimized blade,the airflow mixing downstream of the blade's trailing edge weakens,which is more observable at the middle and tip sections,and hence the entropy increase,which is in the middle and tip region of the S3 cross-section downstream of the optimized blade,significantly decreases,as shown in Fig.17.

The limiting streamline on the suction surface before and after optimization are compared in Fig. 18. As seen in Figs.16(b),(c)and 18,the shockwave position is pushed back at the middle and tip sections of the optimized blade,which reduces the separation area of airflow after the shockwave,and thus the separation line of the optimized blade is pushed back significantly.And,with increasing radial height,this phenomenon becomes increasingly observable.

4.4.Stage 35 optimization verification

The new optimization method is verified by a single-stage transonic axial compressor Stage 35.The Bezier surface parametric control method is applied to the suction surface of the rotor and stator blade respectively,with a total of 18 design variables.The design space of the optimization variables of the suction surface of rotor and stator blade is set to[-1.0,0.5],and the negative values represent the expansion direction of the blades'suction surfaces,while positive values represent the adduction direction.

The rotor blade tip clearance is set to 0.4 mm,and the stator blade root clearance is 0.4 mm.The rotor blade grid number is set to 680000,and the stator blade grid number 750000,which meets the requirement of grid quality.The flow field will be calculated with 1.02 million grids for rotor blade and 1.1 million grids for stator blade after the optimal blade is obtained.We set the scale of the bee colony in IABC algorithm to 80,maximum number of iteration to 3,and maximum number of exploitation to 3.

With the same computer configuration as Rotor 37 verification,in this case of parallel computation using 8 CPUs,each generation of the bee colony optimization takes 24 h and a total of three optimization iterations are performed.Hence,the total runtime is 72 h,within the acceptable range of engineering.The history of optimization cycle is shown in Fig.19.

Fig.15 Comparison of relative Mach distributions at hub,middle and tip sections before and after optimization of Rotor 37.

Fig.16 Comparison of static pressure distributions at hub,middle and tip sections before and after optimization of Rotor 37.

4.4.1.Comparison and analysis of optimization results

Table 3 compares the performance of design points at the design speed before and after optimization.After optimization,the flow of the design points increased by 0.45%,the adiabatic efficiency increased by 0.54%,whilst the surge margin expanded by 1.11%.

Fig.17 Comparison of entropy distributions at S3 section before and after optimization of Rotor 37.

Fig.18 Comparison of limiting streamlines on suction surface before and after optimization of Rotor 37.

Fig.19 History of optimization cycle of Stage 35.

The aerodynamic performance of the Stage 35 compressor is compared at the design speed before and after optimization,as shown in Fig.20.In Fig.20(a),the optimized blade flow is,on average,0.5%greater than that of the original blade within the full operating range,with the same total pressure ratio condition.In Fig.20(b),the optimized blade adiabatic efficiency is 0.3%higher than the average of the original one in the full flow range.Compared with the original blade,the optimized blade flows more smoothly under high back pressure,which improves the surge margin of the optimized blade.

From the comparison of the Stage 35 dynamic and static blade three-dimensional geometry before and after optimization shown in Fig.21,it can be seen that optimized Rotor 35 blade is adducted in the range of 15%or more in the radial direction.The adduct region gradually increases from a range of 35%chord length from the leading edge and 15%of the radial height,to a range of 80%chord length from the leading edge at the blade tip(Fig.21(a)).The Stator 35 suction surface is expanded as a whole(Fig.21(b)).

Fig.22 shows the change of the suction surface of the blade before and after optimization of the dynamic and stationary blades.Fig.23 shows the comparison of the changes of the hub,middle and tip of the Stage 35 before and after optimization(magnification of the geometric change by 10 times).As can be seen from Figs.22(a)and 23(c),the maximum extent of adduction of Rotor 35 blade locates at 35%from the leading edge at the tip,and the maximum extent of external expansion locates at 30%from the trailing edge at the middle of the blade.As can be seen from Fig.22(b),the maximum area of the expansion is 45%of the chord length from the leading edge at the tip,which is verified in Fig.23(d)-(f).

4.4.2.Comparison and analysis of flow field before and after optimization

As shown from Fig.24,the positive incidence angle of incoming airflow of the optimized rotor and stator blade is smaller than that of the original one from hub to tip,and thus theaccelerated distance,where air flows around the leading edge of the optimized rotor and stator blade,shortens.Therefore the intensity of oblique shockwave weakens there,which is verified in Fig.25(b)and(c),reducing the loss of shockwave and that of interaction between the shockwave and the boundary layer.According to Figs.23(b),(c)and 26(c)-(f),the profiles of the middle and tip sections of the rotor blade are Sshaped,and the profile of the first half is adducted,which slows the airflow acceleration at the suction side of the middle and tip sections of the optimized rotor blade,and thus the intensity of channel shockwave weakens,leading to the reduction of the loss of shockwave,whilst the shockwave position is pushed back.Consequently,the air separation area near the trailing edge due to an adverse pressure gradient decreases,reducing the loss of air separation,and thus the low-entropy area downstream of the trailing edge,shown in Fig.27(a)and(b),increases,decreasing the loss of mixing airflow there.

Table 3 Comparison of Stage 35 performance at design point and design speed before and after optimization.

Fig.20 Comparison of aerodynamic performance before and after optimization of Stage 35.

Fig.21 Comparison of geometry of suction surface on Rotor 35 and Stator 35 blades before and after optimization of Stage 35.

Fig.22 Distribution of changes on suction surfaces of Rotor 35 and Stator 35 blades.

Fig.23 Comparison of geometry of hub,middle and tip sections of Rotor 35 and Stator 35 blades before and after optimization(multiply magnitude of change by 10).

Fig.24 Comparison of radial distribution of incoming airflow incidence angles of Rotor 35 and Stator 35 blades before and after optimization.

Fig.25 Comparison of static pressure at hub,middle and tip sections of blades before and after optimization of Stage 35.

Fig.26 Comparison of relative Mach number at hub,middle and tip sections of blades before and after optimization of Stage 35.

Fig.28 shows the comparison of the limiting streamlines on the suction surface of the rotor and stator blades before and after optimization of Stage 35.It can be seen from Fig.28 that the separation forms of the suction surface are all open separation before and after the optimization.Since the suction sides of the middle and tip sections of the optimized rotor blade are S-shaped and the profile near the leading edge is adducted,the shockwave position pushes back,and then the separation position of the boundary layer caused by the inverse pressure gradient pushes back, thus pushing back the separation line of the optimized rotor blade.Therefore the separation loss is reduced.The separation on the suction surfaces of stator blade is angular separation both before and after optimization.Since the suction side of the tip area of the optimized stator blade has the greatest expansion,a pressure gradient that points to the middle of the blade is formed.Thus,the radial migration of airflow in the region from the tip to the middle of the blade enhances,slightly increasing the open separation area, leading to the slight increase of the flow loss.However,compared with the loss reduction of shockwave,and that of the interaction between the shockwave and boundary layer,and that of the separation of the boundary layer,the increased loss of airflow,caused by the increased open separation area,is small.

5.Conclusions

Combining Bezier surface parametric control method and IABC algorithm with verified CFD numerical simulation,a new engineering global optimization method is constructed,and a single-row transonic rotor(Rotor 37)and a singlestage transonic axial compressor(Stage 35)are used to verify this method.The conclusions are as follows:

(1)The Bezier surface parametric control method is successfully applied to the optimization of Rotor 37 and Stage 35,and only 9 optimization variables are used to control each blade,significantly reducing the number of optimization variables,which saves a lot of calculation costs.Moreover,the fine adjustment of suction surface of blades is well realized,with integrity and smoothness.

(2)In the IABC algorithm,the employed bees and the onlooker bees utilize the global and local optimal information when exploring the food source,respectively,replacing the original way of randomly acquiring the information,and the overall exploration information of the food source can be better utilized.Compared with the general GA and ABC algorithms,the capability of global optimization and convergence speed of the IABC algorithm are significantly improved.

(3)The surface parametric control and global optimization methods are successfully applied to the aerodynamic optimization of the transonic axial flow compressor Rotor 37 and Stage 35,and the optimized results are obtained in shorter time(24 h and 72 h,respectively),which is easy to accept in engineering.After optimization,under the constraint that the flow rate is not reduced and the pressure ratio is basically unchanged,the adiabatic efficiency of Rotor 37 at design point increases by 1.02%,whilst the surge margin increases by 0.84%;the adiabatic efficiency of Stage 35 at design point increases by 0.54%,whilst the surge margin increases by 1.11%.

Fig.27 Comparison of entropy distribution at S3 section downstream of rotor and stator blades before and after optimization of Stage 35.

Fig.28 Comparison of limiting streamline on suction surface of rotor and stator blade before and after optimization of Stage 35.

(4)The new optimization tool composed of the Bezier surface parametric control and global optimization methods presented in this paper has the advantages of saving a lot of calculation cost and fast global optimization for aerodynamic optimization of 3D blade of axial compressor,which leads to broad application prospects for optimization of compressor blades.

Acknowledgements

This study was supported by the National Natural Science Foundation of China(No.51576007)and Civil Aircraft Special Research of China(No.MJZ-016-D-30).