Rotor Cascade Shape Optimization With Unsteady Passing Wakes Using Implicit Dual-Time Stepping and a Genetic Algorithm

An axial turbine rotor cascade-shape optimization with unsteady passing wakes was performed to obtain an improved aerodynamic performance using an unsteady flow, Reynolds-averaged Navier-Stokes equations solver that was based on explicit, finite difference; Runge-Kutta multistage time marching; and the diagonalized alternating direction implicit scheme. The code utilized Baldwin-Lomax algebraic and κ-e turbulence modeling. The full approximation storage multigrid method and preconditioning were implemented as iterative convergence-acceleration techniques. An implicit dual-time stepping method was incorporated in order to simulate the unsteady flow fields. The objective function was defined as minimization of total pressure loss and maximization of lift, while the mass flow rate was fixed during the optimization. The design variables were several geometric parameters characterizing airfoil leading edge, camber, stagger angle, and inter-row spacing. The genetic algorithm was used as an optimizer, and the penalty method was introduced for combining the constraints with the objective function. Each individual’s objective function was computed simultaneously by using a 32-processor distributedmemory computer. The optimization results indicated that

The major cause of flow unsteadiness in the interaction between stator and rotor is the combination of circumferencedirected, nonuniform flow caused by the wake of the stator and the relative motion of the rotor. This study focused on the possibilities of improving the aerodynamic performance of a rotor cascade subjected to unsteady flow due to the wakes of the stator cascade located upstream. The study of rotor cascade shape optimization consisted of the following three procedures: stator and rotor cascade analysis, airfoil parameterization, and airfoil shape optimization. To numerically simulate the unsteady flow for a rotor cascade, the unsteady flow-field analysis code was developed (Lee, 2000). The Navier-Stokes solver was based on the Runge-Kutta explicit time marching (Swanson and Turkel, 1992) and diagonalized alternating direction implicit (DADI) schemes (Pulliam and Steger, 1980). A DADI scheme was used for the results in this paper because it was much faster. Local time stepping and implicit residual smoothing (Turkel et al., 1996) were added to the Runge-Kutta explicit time marching scheme to accelerate the convergence rate. To avoid the expense of solving a block matrix equation, the block penta-matrix was diagonalized and reduced to a scalar 354 E. S. LEE ET AL.  From Binder et al., 1984;Fottner, 1990. penta-diagonal matrix equation. For robustness, multiblock capability was incorporated. The algebraic (Baldwin-Lomax) and two-equation (κ-ε) models were tested for turbulence closure (Sahu and Danberg, 1986), and κ-ε was preferred. The implicit dual time stepping method (Arnone et al., 1993) was selected as an unsteady formulation and implemented to both the Runge-Kutta and DADI schemes. The validation of the unsteady analysis codes was performed by four test cases and the calculation results were successfully compared with analytical solutions, experiments, and other numerical results (Lee, 2000). A multigrid method was used for the stator cascade flow calculations only.
To analyze the rotor cascade aerodynamics with passing wake effect, it is necessary to obtain the wake information as an inlet boundary condition for the unsteady rotor analysis. An existing single-stage Deutsche Forschungs fuer Luft und Raumfahrt (DFVLR) turbine cascade (Table 1) was simulated numerically and compared with experimental data.
The steady flow of the stator cascade of airfoils was simulated to obtain the flow information of its wakes (Figure 1). An Htype grid was used, and 32 grid points were placed between the inlet and the leading edge, 129 points on the airfoil surface, and 96 points on the wake line. Nonreflecting boundary conditions were used at the exit, thus allowing for nonuniform exit flow. The exit Mach number was set to 0.76 and the Reynolds number was 1.2 million, based on stagnation properties.
To optimize the rotor cascade airfoil, a parameterization method of airfoil shapes was necessary. Pritchard's (Pritchard, 1985) formulation for geometric shape parameterization was selected because of its simplicity and robustness. It allows for variation in leading edge radius, leading wedge angle, blade inlet angle, and tangential chord length ( Figure 2). These four airfoil parameters were used as design variables during the optimization.

FIGURE 1
Calculated isentropic Mach number distribution on the cascade surface and experimental data from Fottner (1990).

ROTOR CASCADE WITH PASSING WAKE ANALYSIS
At 2.85 mm axially downstream from the stator trailing edge, the flow information (total pressure, total temperature, and flow angle) was interpolated from the steady stator flow-field solution (Lee, 2000). This interpolated flow information was then used as the inlet absolute frame total conditions for the passing wake simulation of the rotor's flow field. The relative translational speed of the rotor cascade was added to the vertical component of the rotor cascade inlet flow velocity. The definitions of lift coefficient, Cl, and total pressure loss, dp t , were The density and velocity at the rotor cascade exit were used as the reference conditions, and S r was the chord length of the rotor cascade. The reference total pressure (in absolute frame) was the stagnation pressure at the inlet of the stator cascade. Average total pressure (in relative frame),p t , was calculated by integrating the total pressure along the y direction at the rotor trailing edge, where y p means the pitch length of the rotor cascade. The total (relative) pressure loss was calculated at 10% of chord length away from the trailing edge in the streamwise direction. The trailing edge loss was included, but not as a design parameter. In the H-grid used for the rotor cascade, 25 grid points were placed between the inlet boundary and the rotor cascade leading edge; 61 points were used on the rotor airfoil surface; and 43 points were used on the rotor wake. The inlet boundary grid points were equally distributed in the pitchwise direction. Dual time-stepping (Arnone et al., 1993) and preconditioning (Turkel et al., 1996) were used. One period of the passing wake was composed of 10 real-time intervals. Each real-time interval required about 100 pseudo-time subiterations in order to reduce the average residual by four orders of magnitude. After 15 periods, a periodic flow condition was reached. Figure 3 shows the

FIGURE 2
Four geometric parameters that represented optimization design variables in the DFVLR rotor cascade. lift and total pressure loss variation with time for three values of axial gap (Table 2). Figure 4 illustrates the perturbation vector plot and entropy contours at four instances, where T is a time period during which a wake passes one pitch length. The perturbation velocity vectors can be calculated by subtracting the time-averaged velocity from the instant velocity. The passing wakes can be seen at the inlet boundary due to the total pressure variation across the wake. These wakes start bowing as they approach the leading edge of the rotor cascade because of the airfoil circulation. The wake is chopped into two parts by the airfoil leading edge. The chopped wake experiences prolongation while it passes along the suction side. These bowing, chopping, and prolongation characteristics of a passing wake are the typical occurrences while it passes along the rotor cascade (Ho and Lakshminarayana, 1994).
The negative jet effects can be found in the wake region due to the defect of velocity magnitude. This negative jet becomes a sequence of opposite directional vortexes along the suction side in Figure 4. These vortexes are distorted and stretched because

FIGURE 3
Variations in the coefficient of lift and total pressure loss with time in the nonoptimized DFVLR stator-rotor cascade.
of the velocity difference between the pressure and the suction sides.

ROTOR CASCADE SHAPE OPTIMIZATION
To improve the aerodynamic performance of a rotor cascade, three required flow conditions were considered in this study. First, the optimized airfoil should satisfy the mass flow rate, which is specified by the user. The mass flow rate can be controlled by a specified backpressure ratio. Second, the rotor cascade should provide as a high lift force as possible. Finally, the total pressure loss between the inlet and exit planes should be minimized in order to maximize rotor efficiency. However, these three parameters are interdependent and influenced by each other.
In this study, only four geometric design variables were used to represent the airfoil cascade geometry (see Figure 2): the airfoil inlet angle (the higher lift coefficient can be obtained by controlling this angle); the inlet wedge angle (the airfoil thickness and slope of the leading edge region can be controlled by this parameter); the leading edge radius (the location of the stagnation point and the pressure peak at the leading edge can be controlled); and the tangential chord length (thus controlling the pitch length and the stagger angle). The trailing edge radius was kept fixed.
Optimizations were conducted for two different axial distances between the stator and the rotor cascades (gap = 54 mm and gap = 27 mm). Two different objective function formulations and constraints were applied for each axial distance.
In this study, a parallelized version of a standard genetic algorithm (GA) optimizer was developed (Dulikravich et al., 1999) and used on a distributed memory machine for aerodynamic shape optimization (Dennis et al., 1999). If a parallel GA is used and the same number of processors is selected as the population size, the objective function for each member of the population is calculated simultaneously. Then, CPU time of parallel GA = (population size) × (number of generations) ×(CPU time of one CFD run)/(number of processors) All ranges of the optimization parameters (Table 3) were normalized between 0.0 and 1.0. The airfoil exit angle and the axial chord length were fixed during the optimization.

OPTIMIZATION TEST CASES
Most optimization algorithms introduce the constraints as a part of the objective function. If some of the individual designs violate the constraints, penalties are added to the objective function value. Therefore, only the designs that satisfy the constraints can survive during the optimization. We chose to formulate the fitness function, f , as follows (Lee, 2000): Here, β cl and β pt are the positive parameters specified by the user before optimization so that β cl + β pt = 1.0. Here, Cl init-ave represents an average coefficient of lifts for the first generation of the design optimization and p t init-ave is an averaged initial total pressure loss, while p t is defined as The four cases of optimization that were tested are shown in Table 4. Notice that βṁ is much greater than β cl or β pt and represents a user-specified penalty weight, whereasṁ spec is the specified mass flow rate and is an additional penalty for violating the constraints. If an individual design violates the constraints, the specified value is added in its fitness function. If an individual satisfies the constraints, then is zero. Thus, case 1 and case 3 minimize the total pressure loss while the lift cannot be less than a specified value and the mass flow rate is fixed. Case 2 and case 4 maximize the lift, but the total pressure loss cannot be greater than a specified value and the mass flow rate is fixed ( Table 5). The specified values of lift,

FIGURE 5
Convergence history of GA optimization in case 1, axial gap = 54 mm.
total pressure loss, and mass flow rate were obtained from the calculation results of the original DFVLR single-stage turbine.
In the genetic algorithm (Dulikravich et al., 1999), the population size was selected to be 15. A distributed-memory parallel computer with 16 processors (Pentium 400 MHz) was used for optimization (Dennis et al., 1999). One processor was assigned

FIGURE 6
Original and optimized DFVLR rotor linear-cascade airfoil shapes for two axial gap values.
to the main genetic algorithm, and the other 15 processors were used to calculate the objective function of each individual. The maximum number of generations was set to 20. Mutation probability was set to 0.01 and nine bits were used in a string for the variation of design variables. Four design variables were normalized to have a value between 0.0 and 1.0. During the optimization, the grid was automatically generated by the transfinite interpolation method and elliptic partial differential equations. Figure 5 shows the convergence history of the maximum fitness function value. During the first five generations, the maximum fitness function was still negative because all individual designs violated the constraints. Starting with the sixth generation, the fitness function became positive, which means that some of the individual designs began to satisfy the constraints. Figure 6 illustrates the optimized airfoil geometries for each of the four cases. Generally, the inlet blade angle and curvature were increased as in a high-pressure turbine case (Madavan et al., 1999), while the tangential chord length was decreased when compared to the original DFVLR cascade. The effective flow angle of attack was thus increased due to the increase in the inlet blade angle. Figure 7 depicts the instantaneous entropy contours and perturbation velocity vectors for four instants of time during a wake-passing period, T . The negative jet, which Hodson (1998) pointed out, can be found at the inlet region corresponding to the wake induced by the stator. This wake shape changes to a bow
as it approaches the leading edge and results in chopping. The chopped wake on the suction side experiences stretching while it passes along the suction side. This negative jet becomes the two vortexes, which have opposite directions, on the suction side.

FIGURE 8
Variation of lift and total pressure loss with time in an optimized DFVLR rotor linear cascade.
The evolution of the lift coefficient and the total pressure loss in the four optimized cascades are shown in Figure 8. The average total pressure loss in the optimized airfoil is less than in the DFVLR cascade, as expected. Due to the constraints, the  Lee, 2000. lift in case 1 is greater than in the original DFVLR cascade. The average total pressure loss in case 2 is almost the same as in the DFVLR cascade, which is specified. However, the average lift is higher than in the original DFVLR cascade. A comparison of case 1, case 2, and the original DFVLR cascade at gap = 54 mm suggests that the three airfoil cascades have approximately the same total pressure loss (Table 6). However, the average lift in case 2 is about 3.7% higher than the lift in the original DFVLR cascade.
The average total pressure loss in case 3 is less than in the DFVLR rotor cascade. However, the average lift is almost the same as in the DFVLR cascade, which is a constraint. Thus, a decrease of axial gap causes stronger unsteady perturbation of the flow. This strong perturbation results in the increase in the amplitude of lift and in the total pressure variation with time.
Optimization performed with the axial gap reduced in half (27 mm) resulted in the geometry of case 4, which has a smaller leading-edge radius than that of case 3 (see Figure 6). The average total pressure loss in case 4 is approximately the same as that in the DFVLR cascade, which is a constraint. However, the average lift is higher than in the DFVLR cascade (see Figure 8), which is the objective function in this case. Figure 7 depicts the constant entropy contours and perturbation velocity vector field. The lift coefficient in case 4 is about 3.0% higher because the exit flow angle became slightly higher than the airfoil exit angle. The total pressure loss remained the same as in the DFVLR cascade because the trailing edge radius was kept unchanged.
It took about 100 h on a 16 processor Pentium II computer for one optimization case to be performed.

SUMMARY
In this research, the possibility and effectiveness of designing cascade shape using an unsteady Navier-Stokes solver were investigated. The following conclusions can be drawn from this research. The DADI scheme requires about 1/3 of the time needed to converge when compared to the Runge-Kutta scheme. A decrease in the axial gap between the stator and rotor results in an increase in the amplitude of lift and in the variation in pressure loss with time. A decrease in the axial gap also causes an increase in the velocity defect in the incoming wakes and a de-crease in the width of the incoming wakes. This effect causes greater unsteady perturbation of the flow around the rotor.
One of the most important factors in numerical passing-wake simulation is grid clustering. To check the effect of grid density on wake resolution, a doubled grid size was used for the same simulation. The numerical difference in the lift and the total pressure loss between the two grids (129 × 65 and 129 × 129) was within 0.2% and 0.1%, respectively. The optimized cascade increases the lift by about 3.7%, while the total pressure loss and mass flow rate are the same as in the DFVLR cascade. If the numerical error caused by grid clustering is incorporated, it follows that the optimized cascade increases the time-averaged lift by about 3.0%.
The combination of a parallel computer and a genetic algorithm optimization is highly desirable because the individual's objective function at each generation can be calculated simultaneously. In this research, 15 processors were used for the rotor cascade shape optimization and the computing time was reduced to approximately 1/15 of the computing time that would have been required by a single processor. One shape optimization took about 4 days, so if a single processor machine had been used, it would have taken about 2 months for the single optimization.
In parallel genetic optimization, each processor was assigned to run an unsteady flow analysis for a different cascade shape. However, the computing time was different for each case, so the main subroutine had to wait until the slowest flow-field analysis was complete. Thus, a new parallel genetic algorithm is needed that can switch to a new flow-field analysis as soon as the current analysis is completed on each processor. Prichard's (1985) parameterization method of cascade geometry was introduced and proved to be robust and versatile because the user can describe an airfoil analytically by specifying only the parameters that are of fundamental interest. Here, four parameters were of interest and were selected as design variables (the blade inlet angle, the inlet wedge angle, the leading edge radius, and the tangential chord length). However, this geometric representation of the airfoil shape clearly did not offer sufficient geometric flexibility of the airfoil configuration since the optimization could not reduce the total pressure loss significantly while satisfying the constraint that the lift could not be less than the lift of the DFVLR cascade.
Consequently, a more flexible formulation for airfoil geometry is required (Dennis et al., 1999). This would offer more design optimization variables and thus greater opportunity for the optimization algorithm to find and optimize the most influential geometric parameters. Finally, a truly multiobjective constrained optimization method is suggested; it would minimize the total pressure loss and maximize the lift simultaneously (Dennis et al., 2000).