Applications of an Improved Aerodynamic Optimization Method on a Low Reynolds Number Cascade

: The effect of cascade aerodynamic optimization on turbomachinery design is very significant. However, for most traditional cascade optimization methods, aerodynamic parameters are considered as boundary conditions and rarely directly used as the optimization variables to realize optimization. Given this problem, this paper proposes an improved cascade aerodynamic optimization method in which an incidence angle and nine geometric parameters are used to parameterize the cascade and one modified optimization algorithm is adopted to find the cascade with the optimal aerodynamic performance. The improved parameterization approach is based on the N on -Uniform Rational B-Splines (NURBS) method, the camber line superposing thickness distribution molding (CLSTDM) method, and the plane cascade design method. To rapidly and effectively find the cascade with the largest average lift-drag ratio within a certain range of incidence angles, modified particle swarm optimization combined with the modified very fast simulated annealing algorithm (PSO-MVFSA) is adopted. To verify the feasibility of the method, a cascade with NACA4412 and a practical cascade are optimized. It is found that the average lift-drag ratios of two optimal performance cascades are respectively increased by 13.38% and 15.21% in comparison to those of two original cascades. Meanwhile, through optimizing the practical cascade of the Blade D500, under different volume flow rates, the pressure coefficient of the optimized cascade is increased by an average of more than 6.12% compared to that of the prototype, and the average efficiency is increased by 11.15%. Therefore, this improved aerodynamic optimization method is reliable and feasible for the performance improvement of cascades with a low Reynolds number.


Introduction
Blade design is of great importance to the efficiency and properties of turbomachinery. Achieving the aerodynamic design of a blade is a very complex and arduous task due to complicated flow phenomena and the interactions among various parameters [1]. A small geometric change of one blade can lead to a deterioration of the aerodynamic performance of the whole machine [2]. In general, the cascade design method is one of the most popular design methods employed for turbomachinery [3][4]. In this method, the blade is stacked by the sections on the different radiuses and shown in Figure 1a. Additionally, the section is projected onto the plane to form the cascade, as shown in Figure 1b. Therefore, the blade performance is affected by the cascade design. With the rapid development of the computer technique, more and more parameterization methods and Blade optimization algorithms have been proposed and used in the design process of a cascade. This means that the time required to design a new cascade is becoming shorter and the aerodynamic performance of the new cascade is becoming better.
(a) An illustration of blades (b) An illustration of a cascade A good parameterization method can not only use fewer design variables to describe an airfoil accurately, but also rapidly re-construct one airfoil in the optimization process [5,6]. The parameterization methods are usually divided into two categories: The constructive method and the deformative method. Each method has been continuously improved by many researchers [7]. The deformative parameterization method is the simpler of the two methods. In this method, a standard airfoil is deformed to generate one new airfoil, in order to satisfy a certain condition. The Hicks-Henne function [8], radial basis function [9], Bezier function [10], B-Spline function [11], and Non-Uniform Rational B-Splines (NURBS) function [12][13][14] are usually used as deformative functions to generate a new airfoil based on the original airfoil. In particular, NURBS [12] is the most popular function due to its ability of local control and its conics description over the curve. However, these deformative functions are used to generate one new airfoil based on the point coordinates of an airfoil. To relate the airfoil shape to the airfoil geometric feature parameters, the camber line superposing thickness distribution molding (CLSTDM) method [15] was proposed. In this method, several airfoil geometric parameters are used to parameterize the half-thickness distribution curve and the mean camber curve through two polynomials. Then, these two curves are coupled to form a whole airfoil. The feasibility of this method has been proved by several works [15][16][17]. In the CLSTDM method [15,16], the blade contours described by many coordinate points can be transformed into functions controlled by several parameterized variables. Moreover, it is convenient for designers to use this method to parameterize one blade based on their experience. Nowadays, it is widely used in the optimization of turbomachinery. However, the aerodynamic parameters are not considered in the parameterization and still act as the boundary condition of flow field computing.
To rapidly find the airfoil with the optimal aerodynamic property, many intelligent optimization technologies have been used in the process. Among all of the published optimization algorithms, the genetic algorithm (GA) [18][19][20][21] and the simulated annealing (SA) algorithm [10,11,22], as two traditional intelligent optimization algorithms, have been widely used in airfoil optimization. They aim to find the airfoil with the optimal performance precisely. However, much time is required to complete the search, which leads to a poor computational efficiency [23,24]. A new intelligent algorithm, known as particle swarm optimization (PSO) and proposed by Kenney and Eberhart [25], can be used to solve this problem. PSO is a population-based, self-adaptive searching optimization method. The principle of this method is based on animal social behaviors, such as birds' migration. However, it is easy to become trapped into a local extreme value or converged to precocity by the standard PSO. Therefore, researchers began looking for improved methods to solve this problem. Shi [26] used a linearly decreasing inertia weight to balance the global and local searching character. However, the local searching capability of this method was weak. Simultaneously, it is hard to predict the maximum iteration number. Clerc [27] set a constriction factor determined by two learning factors to cancel the boundary limits of velocity, and to balance the global and local searching capability. Hu [28] adopted a stochastic inertia weight to replace the linearly decreasing inertia weight. This method could accelerate the convergence velocity and avoid being trapped into a local best solution. However, these methods still have the risk of local convergence.
Most approaches to parameterization of the cascade have been used with only the help of geometric feature parameters, and the aerodynamic parameters were not referred to. The efficiency of optimization was thus negatively affected. Therefore, in this paper, considering the cascade aerodynamic characters, one aerodynamic parameterization approach for a low Reynolds number cascade is proposed. To find the airfoil with the optimal performance during a certain range of incidence angles, a modified PSO-MVFSA algorithm is studied. Furthermore, two cases, such as the cascade with NACA4412 and the blade of FAN D500, are selected to verify the feasibility of the improved parameterization and optimization method.

CLSTDM-NURBS Method
In this paper, for the CLSTDM method, the pressure side and suction side are obtained through the camber curve superposing the half-thickness distribution curve. This method is used to describe an airfoil and is shown in Figure 2. Due to its good ability to design a complex geometry, the twoorder NURBS function defined by Equation (1) is used to describe the camber curve and the halfthickness distribution curve.
where u is the knot vector, n is the order of NURBS ( =2 n ), i is the mark of the control points and i Q is the control point. To solve the two-order NURBS function, the De Boor algorithm [29], which provides a fast and numerically stable way of finding a point on a B-spline curve with the given u in the domain, is adopted and programed.
For camber parameterization, the camber curve is divided into two two-order NURBS curves (solid line), and these two NURBS curves are respectively controlled by several geometric control points ( 0 1 2 , , Q Q Q 3 4 , , Q Q ) shown in Figure 3.
For parameterization of the half-thickness distribution curve, the curve is divided into three parts, including the leading edge half-thickness, the middle half-thickness, and the trailing edge halfthickness, as shown in Figure 4. The leading edge (LE) part is described by one two-order NURBS curve, the middle part by two two-order curves, and the trailing edge (TE) part by one two-order NURBS curve. Four NURBS curves are controlled by nine geometric control points ( 0 1 2 3 4 5 6 7 8 , , , , , , , , P P P P P P P P P ), which are derived by seven geometric feature parameters of the airfoil, such as the leading edge radius 1 R , the trailing edge radius 2 R , the chord line =1.0 L , the thickness gradient angles 1 2 ,   , and the coordinate of the maximum half-thickness point ( , ) x y T T . As is the case for the camber curve, the half-thickness distribution curve also requires continuity. Therefore, at the three geometric control points 2 4 6 , , P P P , it is necessary to maintain collinearity for the relative lines.   , 0

Improved Aerodynamic Parameterization
Unlike conventional parameterization approaches in which the airfoil is only parameterized with the geometric feature parameters [30][31][32], an improved aerodynamic parametrization approach combining the plane cascade design method and the CLSTDM-NURBS is proposed, in which the incidence angle, i , and nine geometric parameters are used as control variables.
For one specific cascade, it can be assumed for the airflow condition that the inlet flow angle 1  is constant and along the tangential line. In this case, a change of the incidence angle can cause a change of the geometric inlet angle of the cascade, and the variation of the incidence angle i  is equal to that of the geometric inlet angle 1A   . From the definition of the geometric inlet angle, it is clear that the slope of the mean camber curve at the leading point is changed with the changing of the geometric inlet angle. This means that the variation of the incidence angle i  has an indirect influence on the modified value of the leading edge angle ' 1  . The relationship is shown in Equation (4).
Simultaneously, the deviation angle  can also be affected by a change of the incidence angle. In order to determine the deviation angle, one semi-rational formula of the deviation angle was used, based on the plane cascade experiments by Howell [33,34], which was suitable for a low Reynolds number cascade and only applicable in the application conditions ( . The definition formulas of the incidence angle i and the deviation angle  have been proposed in the cascade design process [3,35]. In this paper, one special equation combining a semi-empirical formula [33,34] is shown in Equation (7), where  is equal to 0.5 for the rotor blade. Therefore, if the aerodynamic and geometric parameters i ,  , 1  , 2  , and x B are given, the geometric outlet angle 2 A  can be calculated by Equation (7). Additionally, the revised trailing edge angle ' 2  can be obtained by Equation (8). However, in Equation (7), it is found that the cascade solidity  and the outlet flow angle 2  cannot be determined.
In order to obtain the abovementioned two parameters, the diffusion factor D related to the cascade solidity  , the inlet flow angle 1  , and the outlet flow angle 2  , is introduced as a constraint.
This coefficient can be used to control the aerodynamic load of the cascade. The definition of the diffusion factor is shown in Equation (9).
Therefore, during the process of parameterization, the leading edge angle and trailing edge angle shown in Figure 3 are independent variables and no longer depend on the incidence angle. In this way, the control variables of cascade parameterization are determined, such as

Modified PSO-MVFSA
With the extensive use of standard particle swarm optimization (PSO), some drawbacks have been found, such as the local extreme minimum and the precocity. To solve these problems, many works have been published [26][27][28]. Some only adjusted the change of inertia weight to avoid being trapped into a local optimal solution. Moreover, some used two learning factors to balance the local searching and global searching. In this paper, the control variables of cascade parameterization are grouped as a particle. Additionally, the best particle position corresponding to the optimal cascade is obtained by the modified PSO, as shown in Equation (10) [27].
where V is the particle velocity, subscript i , j is the particle sequence,  is the stochastic inertia weight [28],  is the constriction factor [27], 1 c , 2 c represents two learning factors, 1  , 2  represents the random number uniformly distributed in   i j P is the best position in its flight history, , g j P is the best position in the particle swarms, max  is the maximum stochastic inertia weight, min  is the minimum stochastic inertia weight,  is the variance of the stochastic inertia weight, and 3  , 4  is the random number of the standard normal distribution.
However, it is very difficult to judge whether the results are the local optimal solutions or the global solutions when only using the modified PSO. Due to the probability of the simulated annealing (SA) algorithm jumping out of the local value, it can be used to solve the problem. Furthermore, considering that the low searching velocity of the standard SA algorithm can lead to a greater time consumption, modified very fast simulated annealing (MVFSA) [36] is adopted to search for the global optimal minimum. Since the Cauchy distribution depending on temperature was better than the Gaussian distribution, the coefficient of variable perturbation ij y based on the Cauchy distribution has been redefined by Equation (11). In order to improve the efficiency further, the random probability of the acceptance, r P , has been redefined based on the Boltzmann-Gibbs distribution, as shown in Equation (12).
where u is a random number ranging from 0 to 1 , max T is the initial simulated high temperature, is adopted to evaluate the current energy of the particle.

Verification of Modified PSO-MVFSA
Due to its local optimums and premature results, the Rastrigin function is usually used as a test function. The definition of the function is shown in Equation (13).
In theory, when the independent variable j x is equal to 0 , the global minimum optimum of the function is ( ) 0 j f x  . In Figure 5, three different solutions of the function are respectively shown.
When the independent variable range is   20, 20  , the function looks smooth and monotonic, as shown in Figure 5a. With the increase of the resolution, many small peaks and valleys are observed in Figure 5b,c. This means that the function has many widespread local minima. To prove the feasibility and accuracy of the modified PSO-MVFSA, two other PSO-based optimization algorithms, as two comparison algorithms, are adopted to search for the optimal minimum of the Rastrigin function. As the function is two-dimensional, each particle consists of two variables for three kinds of PSO algorithms. The initial particle number, the total iteration number of PSO, and the total iteration number of MVFSA are set to 30, 80, and 40, respectively. Some other coefficients are listed in Table 1, such as the positive value  , the Markov chain length J , and so on. By means of these searching optimization methods, the optimal results are obtained and shown in Table 2. It can be seen that the optimal function value obtained by PG-PSO is smaller than that of the standard PSO. However, the computing time is longer than that of the standard PSO. For the modified PSO-MVFSA, not only its time consumption, but also its ability for minimum searching, are optimal in comparison with the other two algorithms.

Fitness Function
For improved cascade aerodynamic optimization, nine control variables are used to parameterize the cascade and grouped as a particle. In all flights, the selection of the particle with the best position through the optimum target is very important. The lift-drag ratio of the airfoil is related to its capability regarding the power output and aerodynamic loss [16]. Therefore, the optimization proposition is set up as shown in Equation (14).
where i S is the control variable, L D C C is the lift-drag ratio,   , where i c is the weighting factor and subscript ref presents the original airfoil.

Aerodynamic Optimization Process
During the evaluation of airfoil aerodynamic performance two-dimensional computational fluid dynamics (2D CFD) code which utilized the standard k   turbulence model [37,38] to solve the Reynolds Average Navier-Stokes (RANS ) equations was used, and it was only used to simulate the low Reynolds number airfoils [23,39]. Due to the fact that the lift-drag ratio of one airfoil could be obtained quickly by this code, this 2D CFD code was integrated into MATLAB as a fast solver to evaluate the fitness value of each particle. The whole aerodynamic optimization flowchart is shown in Figure 6. Moreover, the detailed process is introduced in the following steps: (a) Inputting coordinate points of one airfoil, and setting coefficients of the modified PSO-MVFSA; (b) Selecting one angle from the sets of incidence angles and nine geometric variables as an initial particle; (c) Conducting perturbation of nine geometric variables of the initial particle to generate the initial particle swarm with one incidence angle, based on the super Latin square method; constructing airfoil swarm by the improved parameterization method; and evaluating the airfoil fitness value by Equation (15); (d) Finding the best previous position of the particle and the best position of the swarm, recalculating the velocity and the position of each particle by adopting Equation (10), and recalculating the fitness values of the new swarms; is the error function), the data related to the particle are unchanged; if not, replace the data and particles with new data and particles; (f) Repeating steps (d) and (e) until the total iteration of the modified PSO is reached; (h) Putting the swarm particles   i S into the optimization process of MVFSA, conducting perturbation, and re-evaluating the fitness by Equation (10); (i) If 0 E   , the corresponding particle is preserved; if not, the corresponding particle as a basic particle is re-disturbed and re-evaluated; (j) If 0 E   , the re-disturbed particle is retained; if not, the particle is accepted with the acceptance probability equation; (k) Repeating steps (g) (h) (i) until the total iteration of MVFSA is reached; (l) Outputting the particle with the largest function F from the preserved particle swarm at different incidence angles respectively; (m) Selecting the best of the optimal particles by MVFSA as the final optimal particle.

Optimization of a Cascade with an NACA4412 Profile
NACA4412 is adopted as the original basic profile of the cascade and evolved into an advanced airfoil by this improved method. The initial particle number, the total iteration number of the modified PSO, and the total iteration number of MVFSA are set to 30, 100, and 50, respectively. Each cascade is described by the improved aerodynamic parameterization method. The average lift-drag ratio of the airfoil of the cascade is calculated by the flow solver. Moreover, the constant airflow boundary condition is set so that the Reynolds number is as high as  Table 3. In actual airflow conditions, the incidence angle of the cascade is not constant and limited to a certain range. Therefore, the aim of multi-point optimization is not to obtain the best aerodynamic performance under one constant incidence angle, but to reach the best whole aerodynamic performance in the whole working range, referred to in Equation (12). In this paper, the incidence angle is uniformly distributed from 0° to 15° by a step of 1.0°. The aerodynamic performance of each cascade parameterized by one incidence angle and eight geometric variables is calculated by the CFD simulation. Based on Equation (13), the fitness value corresponded to each cascade is figured out. Then, the global minimum value is found by the user of the improved PSO-MVFSA method. NACA4412 airfoil and the optimal airfoil are shown in Figure 7, in which the blue curve is the optimal airfoil and the red curve denotes the original airfoil. From this figure, it can be observed that the suction side and pressure side of the optimized airfoil are changed in terms of geometry in comparison with that of NACA4412. To evaluate the aerodynamic performance of the cascade, the pressure coefficient p C is adopted, which is defined by Equation (16).
where p , 0 p , p  are the current pressure, the stagnation pressure, and the inlet airflow pressure, respectively;  is the density; and u  is the inlet airflow velocity.
The increasing curvature of the suction side close to LE can lead to the intensifying of the velocity increasing and the pressure decreasing. Then, the pressure coefficient of the suction side of the optimized p C near LE is larger than that of the original cascade. However, it is due to the geometric change of the pressure side close to LE that the decreasing of the velocity and the increasing of the pressure are alleviated. These phenomena can be observed in Figure 8. In the figure, the detailed pressure distributions on the suction and pressure side are shown, which respectively correspond to four incidence angle conditions of 0 i   , 3 i   , 6 i   , and 12 i   . It is also clear that the optimized airfoil surface pressure distribution presented by the blue curve is better than that of the original airfoil described by the red curve under each condition. Therefore, the pressure differences between the pressure side and suction side of the optimized airfoil are larger than those of the original airfoil. The cascade performances under the incidence angles ranging from 0° to 15° are shown in Figure 9. The average lift-drag ratio of the optimized cascade is increased by 13.38% in comparison with that of the cascade with NACA4412.
For further analysis, CFD simulation software CFX was used to calculate the cascade performance. The velocity contours of the original and optimized cascades at incidence angles of 0 i   and 12 i   are shown in Figure 10. It is clear that the velocity difference between the suction side and pressure side of the optimized cascade is bigger than that of NACA4412. Therefore, the diffusion factor D of the optimized cascade can be increased. It can also be observed that the airflow separation point on the suction side of the optimized cascade moves backward along the airfoil in comparison with that of NACA4412 and the area of the wake becomes smaller than that of NACA4412. Therefore, it can be considered that the aerodynamic load of the cascade is increased and

Blade D500 Optimization
Blade D500 is used in an axial-flow fan for an evaporator system. It is a kind of low Reynolds number 3D blade. In this work, in order to verify the feasibility of the improved method applied to the practical blade design, the cascade at the 50% radius of Blades D500 was selected and optimized. Two performance parameters, including the pressure coefficient p C defined by Equation (16) and the efficiency  defined by Equation (17), were used to evaluate the aerodynamic performance of Blade D500. In order to obtain the aerodynamic performance parameters, CFX was used to calculate the airflow field of Blade D500.
where v Q is the volume flow rate, P is the static pressure increase, and N the aerodynamic power.

Validation of the CFD Simulation Based on Experiments
The experiments of the axial fan with Blade D500 shown in Figure 11 were conducted in the Key Lab for Power Machinery and Engineering of SJTU. In order to match the practical situation, the rotor was mounted on the guard grill, which was connected with a short bell month at the inlet of the tube. The diameter of the tube was increased to 600 mm to avoid destroying the flow field at the rotor outlet and simulate the real situation in the evaporator system as much as possible. The pressure probes were mounted at points A and B, by which the flow rate and the pressure increase could be calculated. Under eight volume flow rates, the pressure increases were calculated. Under the same conditions as those used in the experiments, CFX was adopted to simulate the aerodynamic performance of the axial fan. The grids consisted of the structured hexahedral meshes shown in Figure 12 that were generated by TurboGrid. The total number of nodes was 1,260,626 for the full channel after studying the grid independence, and the minimum value of the mesh quality was 0.3. Due to the excellent ability to simulate the flow with a fiercely adverse pressure gradient of the standard k   model, it was selected as the turbulent model to simulate the flow field. The mass flow and static pressure were set at the inlet and outlet, respectively. After a series of experiments and CFD simulations of the axial fan with Blade D500, the pressure coefficients were obtained and are shown in Figure 13. Based on Figure 13, it could be found that the pressure coefficients sp C calculated by CFD are slightly higher than those calculated by experiments.
However, the trends of the two curves are very similar. In Figure 14, the efficiency is compared for values obtained by calculations and experiments. It can be observed that the relative error between the two series of average efficiencies is under 8.35%. Considering the abovementioned results, the CFD simulation can be used to feasibly estimate the aerodynamic performance of an axial fan.

Optimization of Blade D500
For aerodynamic optimization of the cascade at the 50% radius of Blade D500, the airflow condition was set as After a series of optimization iterations, the airfoil of the optimal cascade was obtained and is shown in Figure 15, as well as the airfoil of the original cascade. It can be observed that the airfoil of the optimized cascade is different from that of the original cascade. These geometric changes of the optimal cascade can lead to an increase of the curvature of the suction side, while the pressure side is changed a little. Additionally, aerodynamic performance comparisons are shown in Figure 16. It could be found that the pressure differences of the optimized cascade are larger than those of the original airfoil along the chord line within the range of the incidence angle. After a series of simulation calculations, it could be found that the average lift-drag ratio of the optimized cascade is increased by 15.21% in comparison with that of the original cascade. Therefore, it can be considered that the aerodynamic performance of the optimized cascade is better than that of the original one.  To evaluate the feasibility of the improved aerodynamic optimization method, the optimized 3D Blade D500 and the original 3D Blade D500 were simulated by CFD under the same airflow condition. Moreover, the simulation results of two blades were compared. The streamlines of the suction surfaces of the two blades are shown in Figure 17. From the figure, it is clear that the streamline between the middle section and the hub of the original blade is not very good, and a large turbulence loss must have occurred in this region. In contrast, the streamline in the same region of the optimized blade becomes very smooth. The pressure coefficients sp C of the two fans are shown in Figure 18. In this figure, it can be seen that the pressure coefficient of the optimized fan is better than that of the original fan, and the average pressure coefficient is increased by 6.12%. Meanwhile, from Figure 19, it can be observed that the efficiency of the optimized fan is higher than that of the original one, and averagely increased by 11.15%. Therefore, from the above analysis, it can be considered that the aerodynamic performance of optimized Blade D500 is improved and the improved aerodynamic optimization method is feasible.

Conclusions
In this paper, an improved optimization method especially applied to a cascade with a low Reynolds number is proposed. In this method, the incidence angle and eight geometric parameters are used as the control variables altogether. Meanwhile, the modified PSO-MVFSA algorithm is adopted to optimize the objective cascades, and two cascade cases, such as NACA4412 and Blade D500, are selected to testify the method. Some valuable results are as follows: (a) Since the aerodynamic parameter, such as the incidence angle, is considered as one of the control variables, the relationship between the geometry of the airfoil and the aerodynamic performance of the cascade is learnt. Therefore, during the whole optimization process, an improvement of the aerodynamic performance can give rise to a direct modification of the geometry so that the optimization becomes more targeted and more efficient; (b) In this study, particular effort was devoted to designing a fitness function which is suitable for optimizing a cascade with a low Reynolds number. Furthermore, the combination of PSO and MVFSA succeeded in increasing the optimization efficiency and avoided the local optimal to reach a global solution, which was verified by the Rastrigin function and two cascade cases; (c) Based on the analysis of the results from the two cascade cases, such as NACA4412 and Blade D500, it was demonstrated that the average lift-drag coefficient of the optimized cascade was improved, whilst the drag coefficient was kept at a low level. Therefore, it can be considered that the modified PSO-MVFSA can be adopted as an efficient and robust optimizer to solve the problems of multi-variable optimization confronted in cascade design.
Author Contributions: S.Z. and B.Y., investigation, methodology, software, and writing-original draft preparation. H.X. and M.S., validation. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.
Acknowledgments: This research was supported by the National Science and Technology Major Project of China (2017-Ⅱ-0006-0019).

Conflicts of Interest:
The authors declare that they have no conflicts of interest with regards to this work. They declare that they do not have any commercial or associative interests that represent a conflict of interest in connection with the work submitted.