Adjoint aerodynamic optimization of a transonic fan rotor blade with a localized two-level mesh deformation method

Accepted Manuscript Adjoint aerodynamic optimization of a transonic fan rotor blade with a localized two-level mesh deformation method

Xiao Tang, Jiaqi Luo, Feng Liu

PII: DOI: Reference:

S1270-9638(17)30647-8 https://doi.org/10.1016/j.ast.2017.11.015 AESCTE 4290

To appear in:

Aerospace Science and Technology

Received date: Revised date: Accepted date:

12 April 2017 5 August 2017 7 November 2017

Please cite this article in press as: X. Tang et al., Adjoint aerodynamic optimization of a transonic fan rotor blade with a localized two-level mesh deformation method, Aerosp. Sci. Technol. (2017), https://doi.org/10.1016/j.ast.2017.11.015

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Adjoint Aerodynamic Optimization of a Transonic Fan Rotor Blade with a Localized Two-Level Mesh Deformation Method Xiao Tanga,1,∗, Jiaqi Luoa,2 , Feng Liub,3 a Department

of Aeronautics and Astronautics, Peking University, Beijing, China of Mechanical and Aerospace Engineering, University of California, Irvine, United States

b Department

Abstract We present an optimization platform for turbomachinery with complex mesh configuration in a parallel computation environment. A continuous adjoint solver for 3-D viscous internal flow is coded under the same parallel framework as the flow solver. To meet the various permitted extents of reshaping on blade surface and to cut down the computational cost in grid perturbation, a localized two-level mesh deformation method is developed based on Gaussian radial basis function (RBF). This method works efficiently for both the O mesh surrounding the blade and the O-H mesh inside tip gap. In the optimization of the transonic NASA Rotor 67 for high adiabatic efficiency with a mass flow rate constraint, an adjoint sensitivity analysis is conducted. The relations between the design sensitive regions and physical phenomena in internal flow are discussed. Flow fields before and after the adjoint optimization are investigated, including shock system, tip leakage flow, and flow separation. Keywords: Adjoint Method, Aerodynamic Design Optimization, Mesh Deformation, Compressor Rotor, Transonic

Nomenclature

Bs

Source matrix due to rotation

α

RBF bump decaying rate

ξ

Weight coefficient of the mass flow rate constraint

β

Flow turning angle

a

Local sound speed

w

Flow variable vector

BIO

Inlet and outlet boundaries

m ˙

Mass flow rate

cp , c v

Specific heat at constant pressure, volume

η

Adiabatic efficiency

G

Adjoint gradient

Γ

Geometric boundary

h

RBF bump height

u ˆi

Grid velocity components

I

Cost function

κ

RBF bump shape factor, κ = h/λ

Ni

Unit normal vector components

λ

RBF bump half decay radius, λ = ln2/α

p, ρ

Pressure, density

sgen

Total entropy generation

Sij

Projected areas

ui

Velocity components

λi , ki Eigenvalues, eigenvectors of matrix A1 π, θ

Total pressure ratio, total temperature ratio T

Ψ

Adjoint variable vector, Ψ = (ψ1 , ψ2 , ψ3 , ψ4 , ψ5 )

c

Vector comes from the cost function at the outlet

ˆ i Ie As,i , Ai Jacobian matrices, As,i = Ai − u ∗ Corresponding

author Candidate 2 Research Associate 3 Professor, Fellow AIAA 1 Ph.D.

Preprint submitted to Elsevier

1. Introduction Turbomachinery blades are characterized by highly sensitive areas such as the leading edge, suction surface shock region and blade tip, which directly relate to complex flow phenomena like the shock-boundary layer interaction, tip leakage flow, shock-tip leakage vortex interaction. These areas have various geometries as well as levels of design November 13, 2017

sensitivity and permit different extents of reshaping in design optimization. Thus, it is significant to make sure that localized geometry parametrization is used in design and optimization. However, in the widely used shape function perturbation approaches[1, 2, 3, 4, 5], one shape function usually covers a too large area of the blade such that localized sensitivities with high resolution in key design regions are impossible. Discrete sensitivity calculation on each blade surface grid node has also been tried[6], but is very time consuming and requires additional grid smoothing which harms the sensitivity accuracy. Automatic design optimization for turbomachinery components under complex internal flow condition is difficult, especially when large flow separation and strong tip leakage flow are present. The reasons are two-fold. Firstly, CFD codes face difficulties in predicting internal flow separation and shock-vortex interaction. Secondly, under such nonlinear flow conditions, the optimization target functions are probably far from any continuous one. Though adjoint method has been applied to many scientific/engineering problems, including subsonic intake shape optimization[7], multi-row turbomachinery optimization with mixing plane[8], multidisciplinary design optimization of aerodynamics and aero-elastics[9], and optimization in unsteady cascade flow with harmonic balance method[10], its physical contents are far from clear. On the one hand, the interpretations of the adjoint solutions are needed for different kinds of optimization targets which require different kinds of boundary conditions for the adjoint fields. Marta et al.[5] tried to interpret the adjoint variables as the sensitivity of the corresponding conserved flux quantities in the Navier-Stokes equations to metrics of interest. On the other hand, bridges between the adjoint sensitivities and physical phenomena should be constructed, especially in complex flow conditions. In this study, the ability of the adjoint method as an automatic optimization tool under highly nonlinear complex flow condition in turbomachinery is explored. Firstly, the authors present a localized two-level mesh deformation method in a multi-block parallel environment, which not only provides good control for design sensitive regions like blade tip and leading edge but also requires a small amount of computation time. Then based on this method, a high-speed transonic fan rotor blade row with flow separation and strong tip leakage flow near the working point is optimized using a continuous adjoint method. The optimization target is to increase the adiabatic efficiency or equivalently to reduce the entropy generation at a point near peak efficiency. The mass flow rate constraint is introduced by a penalty function approach with a weight coefficient to balance the constraint and optimization target. Firstly we conduct an adjoint sensitivity analysis for the whole blade surface. Physical interpretations are given for design sensitive regions like the shock wave and blade tip. Blade geometry before and after optimization are compared. The influences of blade reshaping upon flow details 2

are also discussed. It is interesting that small modifications on the blade surface geometry at tip considerably vary the tip clearance flow and consequently changes the blockage caused by tip leakage vortex-shock interaction in the flow passage. The weakening of shock system favors increasing adiabatic efficiency through lessening shock loss and flow blockage, as well as suppressing the flow separation at outer span on the suction side. 2. Adjoint Method In 1988, Jameson[11] introduced the adjoint method in airfoil design optimization. This method was successfully expanded to turbomachinery design optimizations[12, 13, 2, 3, 4, 9, 14, 8] in the recent decades. Currently, both the continuous and discrete adjoint methods are used by researchers. In the continuous approach, the adjoint equations are obtained by variation of the flow governing equations and then discretized, while the discrete approach formulates the adjoint system directly from the discretized flow equations[14]. The continuous adjoint method is adopted in this work because it requires less computational resource and the same numerical techniques used in the flow solver can be applied directly to the adjoint solver. 2.1. Adjoint Principle Generally, in aerodynamic design optimization, the cost function I depends on both flow variables w and geometric boundary Γ, i.e., I = I (w, Γ). Variation of the cost function is  T   ∂I ∂I δI = δΓ (1) δw + ∂w ∂Γ In steady state, the flow variables and geometry are connected through the residual formula R (w, Γ) = 0, which by variation gives T    ∂R ∂R δΓ = 0 (2) δw + δR = ∂w ∂Γ Introduce adjoint variable Ψ and rearrange.  T T       ∂I ∂I T ∂R T ∂R −Ψ δΓ δw+ −Ψ δI = ∂w ∂w ∂Γ ∂Γ (3) Then eliminate the variation of flow variables in Eqn. (3) by solving the adjoint equation     ∂I ∂R − Ψ=0 (4) ∂w ∂w such that δI only depends on δΓ, i.e.,     ∂I ∂R − ΨT δΓ δI = ∂Γ ∂Γ

(5)

Thus the gradient with respect to the geometric parameters is     ∂R ∂I δI = − ΨT (6) G= δΓ ∂Γ ∂Γ

The eigenvalues of matrix A1 are λ1 = u1 − a, λ2,3,4 = u1 , λ5 = u1 + a, here u1 and a are the axial velocity and local sound speed. The corresponding eigenvectors are ki , i = 1, 2, 3, 4, 5. On the inlet/outlet boundaries, along each characteristics, the characteristic equation kiT ·dΨ = 0 holds. This leads to

2.2. Adjoint Equations and Boundary Conditions The adjoint system of equations and their corresponding boundary conditions for the inviscid and viscous flow of both stators and rotors have been described in references[2, 3]. Here a brief summary are given together with some improvements on the inlet/outlet boundary conditions. For rotating turbomachinery, the viscous adjoint system is formulated as ∂Ψ 1  −1 T M + Y + BT (7) AT s,i sΨ=0 ∂xi J

kiT · Ψ = constant

(11)

since ki do not depend on adjoint variable Ψ. Equations (10) and (11) together provide the non-reflection inlet/outlet boundary conditions for the adjoint field. A linear system of five equations needs to be solved for Ψ on each cell surface at the inlet/outlet boundaries.

where As,i = Ai − u ˆi Ie , Ai and Ie are the Jacobian matrices and the unit matrix, while u ˆi are the grid velocity T  components. 1/J M−1 Y is the viscous adjoint operator, which comes from a combination of the adjoint parts of the momentum and energy equations. BT s Ψ represents the source terms due to rotating effect. The thermal conductivity and viscosity coefficients are assumed constant in the formulation of the adjoint equations. This simplification is commonly used because, on the one hand, the geometry modification and the flow field variation are small in each design optimization step, on the other hand, since the Navier-Stokes equations and the turbulence model equations are solved in every step, the thermal conductivity and viscosity coefficients are also updated. The frozen eddy viscosity assumption is also used such that no adjoint equation for the turbulence model equations is needed. On the adiabatic wall, the viscous boundary conditions are ∂ψ5 = 0, ψi+1 + u ˆi ψ5 = 0, i = 1, 2, 3 (8) ∂n At the inlet/outlet boundaries, the viscous effects are neglected. The continuous adjoint method requires the integral

 T c − Ni Sij ΨT As,j δwdB (9)

3. Description of Design Optimization NASA Rotor 67 is an axial-flow transonic fan rotor with 22 blades and an aspect ratio of 1.56. At the design point, its tip relative Mach number is 1.38, pressure ratio is 1.63 at a mass flow of 33.25 kg/sec with a rotational speed of 16043 rpm. The running tip clearance is about 0.625% of the leading edge blade height. Strazisar et al.[15] gave detailed laser anemometer measurements. The multi-block grid for Rotor 67 used in this paper is shown in Fig. 1(a). The blade is surrounded by an O mesh, which has been split into three sub-blocks along the spanwise direction. The y + values of the mesh points just above the pressure/suction surfaces of the blade are 1.0. Fig. 1(b) presents the mesh topology at the blade tip. The tip gap is filled with O-H type blocks with 8 cells in the radial direction. The number of grid nodes in one passage is about 430,000. This grid is relatively coarse, especially inside the tip gap. However, for iterative design optimization, this grid offers an acceptable computational cost while capturing all the important flow features, like shock-boundary layer interaction, flow separation, and tip leakage flow. A parallel 3-D Navier-Stokes solver is used to calculate the internal flow field with the classical JST scheme[16] and SST k-ω turbulence model[17]. The choked mass flow rate given by the computations is 34.59kg/sec while the measured is 34.96kg/sec. Figure 2 shows the computed working characteristics in comparison with experimental data. The adjoint system in Eqn. (7) is solved in the same parallel architecture with similar numerical techniques of the CFD solver (multigrid iteration, spatial discretization, Runge-Kutta time marching method with local time step). Artificial dissipation is also needed in solving the adjoint equations. The dissipation coefficients, as in the flow field solution, are based on local static pressure changes. The optimization target is to achieve higher efficiency or equivalently lower losses at a point near peak efficiency with back pressure equal to 1.04 times of the ambient static pressure. Here we take the total entropy generation, which measures aerodynamic losses, as the optimization target.

BIO

to be zero. Here BIO represents the inlet/outlet boundaries. Vector c comes from the cost function defined at the outlet. Ni and Sij are the unit normal vector components and projected areas in the computational domain. δw is the variation of the flow variable vector. The inlet/outlet boundary planes are perpendicular to the axial direction, i.e., N2 , N3 = 0. So, the inlet/outlet boundary conditions only depend on the Jacobian As,1 , which equals to A1 . To ensure zero value of the integral in Eqn. (9), we let  T (10) c − Ni Sij ΨT As,j δw = 0 at inlet/outlet boundaries, which provides the first part of the inlet/outlet boundary conditions. Linear extrapolations of some selected adjoint variables could be used as the rest part of the inlet/outlet boundary conditions (based on the wave propagation directions of the hyperbolic adjoint system). To avoid wave reflections on the inlet/outlet boundaries, Riemann invariants of the adjoint system are adopted. 3

0.95

2.00

0.90

1.80

0.85

1.60

0.80

π − CFD π − Exp η − CFD η − Exp

1.40

1.20

0.93

0.94

Adiabatic Efficiency

Total Pressure Ratio

2.20

0.75

0.95 0.96 0.97 0.98 Relative Massflow

0.99

0.70 1.00

Figure 2: Working characteristics of Rotor 67

condition is

(a) Computation grid (duplicated and rotated) Z Y X

⎡ ⎤ ⎡ ⎤ 0 c1 ⎢Ni Si1 ⎥ ⎢c2 ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ c = ξ ( m ˙ − m ˙ ) c=⎢ 0 ⎢Ni Si2 ⎥ ⎢ 3⎥ ⎣Ni Si3 ⎦ ⎣c4 ⎦ c5 0

(13)

where c is the vector coming from the cost function. ξ is the weight coefficient of the mass flow rate constraint. Ni and Sij are the unit normal vector components and the projected areas in the computational domain, respectively. 4. Mesh Deformation In the optimization cycles, the crucial step is to obtain the gradient sensitivities with respect to the design parameters, usually by adding a small perturbation to the blade shape. When the finite-difference method is used for the gradient calculation, the number of flow solutions and the number of mesh perturbations are both equal to the number of the design parameters. Introducing the adjoint method in the gradient computation significantly reduces the number of flow solutions, yet it does not cut down the mesh perturbation cost because the finite-difference evaluation is still used to obtain the gradient of the mesh deformation to the perturbation of each design parameter. Modern CFD mesh often consists of millions of grid nodes, making it impractical to generate a new mesh from zero whenever a perturbation is added to the shape. Instead of generating a new mesh, an alternative way is mesh deforming, i.e., passing the geometric perturbation to the current mesh and make sure it provides accurate adjoint gradient (if the adjoint method is used) and maintains reliable grid for the next CFD calculation. In a shape optimization, researchers often use smooth shape functions[1] as the geometry perturbation such that the shape remains smooth during optimization process and the number of design parameters can be easily controlled. However, the disadvantage of using shape functions is: the

(b) Grid at tip Figure 1: Computation grid of Rotor 67

To maintain designed thrust, the cost function I is formulated to be a combination of the total entropy generation sgen and half the square of mass flow rate deviation. 1 2 I = sgen + ξ (m ˙ −m ˙ 0) 2

(12)

where the total entropy generation is sgen = BO ρu (Δs) dB in which Δs = cv ln (p/p0 ) − cp ln (ρ/ρ0 ) is the entropy generation per unit mass flow. The weight coefficient ξ = 104 . The benefit of this penalty function approach is that the original constraint optimization problem now becomes a non-constraint one such that no matter how many constraints there are, only one adjoint calculation is needed for the complete gradient. In the adjoint method, different cost functions lead to different forms of (δI)Γ and (δI)w , which result in different adjoint boundary conditions[2, 4]. With the mass flow constraint in Eqn. (12), the corresponding outlet boundary 4

obtained sensitivities are global, thus it is hard to deal with local highly sensitive regions such as the blade tip or the leading edge. On the other hand, sensitivities to the perturbation of every grid node on the design surface may be adopted. This grid node based perturbation gives highly resolved sensitivities but needs much more computational time and requires additional steps to smooth the optimized geometry[6]. We propose in this paper to use the Gaussian radial basis function (RBF) bump to perturb the geometry. Firstly, because of its rapidly decaying property, accurate local sensitivities can be obtained by placing the RBF bumps at every grid node on the blade surface. Secondly, RBF bumps are continuous functions such that no additional smoothing is needed. The decaying rate of RBF can be easily controlled, as shown in subsection 4.1. Based on the RBF, we develop a localized two-level mesh deformation method for perturbing the blade surface as well as passing the surface perturbation to the inner (into the CFD mesh) grid nodes.

expensive to solve in a reasonable amount of time. Though several measures are adopted to reduce the computation cost, the time required in geometry perturbation is still comparable to flow solution. We propose a mesh deformation method which avoids brutal force computation. The first idea is that mesh nodes can be classified into 1st and 2nd levels during perturbation. As illustrated in Fig. 3, the blue dashed lines represent the original mesh and the black solid lines are the perturbed mesh by one RBF located at P (The perturbation effect has been exaggerated). For node P, the grid points on line PQ (radial line perpendicular to the blade surface) are the 1st level nodes while points on lines that lie orthogonally to PQ are the 2nd level nodes, such as the lines SP and TH and many other lines which may or may not intersect with PQ. Notice that the 1st level nodes are affected directly by P’s movement, then they pass the perturbation (with a decaying magnitude) to the 2nd level nodes through the net of grid lines. The 1st and 2nd node levels include all the mesh nodes. There are actually two ways for interpreting the movement of point T due to the change of the boundary point P. In the first way, perturbation goes along the P-H-T path, i.e., from surface node P to 1st level node H and then to 2nd level node T. The second way is the P-S-T path, i.e., from surface node P to surface node S and then to T, the 1st level nodes of point S. From the second interpretation, it is clear once the perturbation at P is spread smoothly to the points on the blade surface, it is sufficient to perturb only 1st level nodes to pass the blade surface movement to the interior points in the mesh blocks. Accurate gradient evaluation and reliable new mesh for next step CFD computation can both be realized. Let’s go through the stated mesh deformation steps again. Firstly, perturb the blade surface nodes. In this step, the movement of each surface node consists of the influence of all the RBFs on the blade surface. So the spreading from P to S (see Fig. 3) is done. Secondly, traverse all the blade surface nodes and perturb their 1st level

4.1. Radial Basis Function The mathematical expression of the Gaussian radial basis function (RBF) bump is f (r) = he−αr

(14)

where h is the bump height (and also the design parameters), α controls the bump decay rate, r is the distance from the bump center. Define half decay radius λ and bump shape factor κ, λ=

h ln2 , κ= α λ

(15)

where κ is kept constant for all RBF bumps while λ is determined locally. 4.2. Localized Two-Level Mesh Deformation Different regions on the blade surface may correspond to different physical phenomena and thus display different adjoint sensitivities and permit different extents of reshaping. Therefore we use localized RBF bumps for geometry perturbation. RBF bumps are placed at every grid node on the blade surface. The half decay radius λ is determined locally. λ ∼ local cell length

Q

T H S

(16)

P

Then α and h are calculated by Eqn. 15. The total perturbation on a blade surface node is the summarized effect of all the RBF bumps. The next step is to pass the perturbation into the whole mesh block. The perturbing magnitude decays in this process to ensure the far field computation boundaries do not change. Walther and Nadarajah[8] used a full perturbation approach, i.e., for every grid point, consider the influences of all RBFs. This method introduces a matrix which is too

Figure 3: Two-level perturbation: ◦ first level nodes of point P

5

section 4

nodes. Table 1 compares the number of algebraic operations of our localized method and the full deformation method in both gradient calculation and the new mesh generation. Assume the number of mesh nodes in each of the three spatial directions is O (n). The numbers of  grid nodes and RBFs on the blade surface are both O n2 . The total number of grid points is O n3 . When calculating the gradient, the first step of mesh deformation uses one single RBF each time and compute its effect on every blade surface node. The corresponding operational count  grid is O n4 , equal to the number of RBFs times the number of surface nodes. The second step is to pass the blade surface movement into the whole mesh blocks. With our localized two-level method, we traverse each surface node and perturb its 1st level nodes with an operational count  of O n5 = number of RBFs × number of surface nodes × number of 1st level points of each surface node. In the full mesh deformation method, the corresponding operational count grows to O n7 because all mesh nodes (the 1st and 2nd levels) are considered to be influenced directly by all the RBFs on the blade surface. After the gradient calculation, a new mesh needs to be generated for the CFD calculation in the next cycle. Now the movement of every blade surface node are affected by all the RBFs on blade surface, giving O n4 operational count. Then, to pass the surface perturbation blocks, the op into mesh erational counts are O n3 and O n5 for the localized two-level mesh deformation method and the full method, respectively. The movements of blade surface nodes are small compared to the mesh blocks. Also, the decaying perturbation magnitude is small compared to the cell size in the middle flow region. We limit grid perturbation inside the O mesh surrounding the blade, i.e., let the decaying perturbation magnitude to be zero at the outer edge of the O mesh. This will not influence the gradient because an exponentially decaying rate is adopted in passing the surface disturbance into the mesh block such that the perturbation magnitude has already become negligible when reaching the O mesh outer edge.

N E A D

Gradient calculation New mesh generation

G

F

S Q

section 3

M L

P section 2

Figure 4: Tip gap O mesh divided into 4 sections

(a) Leading edge

(b) 75% Chord length

(c) Trailing edge Figure 5: Grid perturbation for tip gap mesh blocks

direction and the movements of their outer edges are determined by perturbations of the inner surface of the O mesh surrounding the blade tip. The O-H blocks are deformed as a whole. To apply our deformation method, the i-j plane of the O mesh is divided into 4 sections, as shown in Fig. 4. The deformation starts from the cutting line AB and goes clockwise. In section 1, based on our localized two-level idea, the movement of node A is passed into the O-H mesh blocks along line segment ABCDE, with an exponentially decaying rate. And the perturbing magnitude of point E is forced to be zero. The disturbance of node F spreads along line FGHSK and ends at point K. For section 2, LMGBA and PQHCN are the 1st level perturbing line. The same operation is conducted for section 3 and 4. The long and narrow O-H type blocks inside the tip gap (see Fig. 1(b)) increase the difficulty in three aspects.

Table 1: Time complexity

Blade Surface  O n4 O n4

section 1

B

H K

4.3. Mesh Deformation in Tip Gap The method proposed in subsection 4.2 is also adopted for tip gap mesh blocks. There are 8 cells in the span

Deformation

C

Mesh Blocks Local Full   O n5 O n7 O n3 O n5

6

Firstly, the O-H type blocks are so narrow that the perturbation on one side must be damped efficiently. Otherwise, its influence might penetrate the mesh blocks (along line PQHCN in Fig. 4 for example) and damage the cell layout on the other side. Secondly, perturbing grid nodes along these lines requires to work on two different blocks which might not be on the same CPU such that additional MPI operations are needed. The third thing is that the adjacent nodes between the O and H mesh must be handled carefully to avoid singularity. Figure 5 shows the perturbation of one RBF on tip gap mesh blocks, including three locations: leading edge, 75% chord and trailing edge, with exaggerated perturbation magnitude. Here the blue dashed lines are the original mesh while the black solid lines are perturbed mesh. Based on a local scaling idea, our method achieves robust perturbation for all blade regions and works efficiently in providing accurate gradient and reliable mesh for next optimization cycle. The proposed two-level mesh deformation method is most convenient to apply in structured grids where it is easy to define a grid node as the first level node of a certain blade surface node. Also, the orthogonally aligned i-j-k lines in structured grids help implement the two-level mesh deformation efficiently. The most important thing to be considered when using the two-level mesh deformation method in an unstructured grid is how to define the first level nodes of a certain blade surface node, such that no grid node would be left out. This issue strongly depends on the node connectivity in wall normal direction and might require elaborate algorithms.

The adjoint gradient is positive in the purple regions and negative in the orange. On SS, the blue dashed line approximates the location of passage shock. Before the shock, adjoint gradient components are positively large, from about 30% span to the blade tip. After the shock, the adjoint sensitivities are negatively large, starting around 17% span and end at the blade tip. This says that to reduce shock loss, the SS needs to be thinned before the shock while thickened after it. This effectively creates an inverted S-shaped suction surface to reduce the passage shock into a series of compression or weaker shock waves on the blade surface. Blade to blade relative Mach number contours on 73.4%, 86.0%, 95.8% and 99.3% span are shown in Fig. 7. At the 73.4% span, the relative Mach number behind the bow shock drops slowly from 1.0 to 0.9. In sub-figure (b) at the 86% span, the bow shock at the leading edge becomes stronger. After it, a passage shock comes into being. This two-shock system has been reported by Strazisar[18] and is typical in transonic rotors. Sub-figure (c) gives the Mach number contours at the 95.8% span, where the shock system is further strengthened such that the bow shock and passage shock stand very close to each other on SS. The last plot in Fig. 7 presents the relative Mach number contours at the 99.3% span where the effect of the tip leakage vortex (TLV) is clear. Pandya and Lakshminarayana[19] and Dawes[20] pointed out that the main part of the tip clearance flow initiates close to the leading edge, rolls up in the flow passage to form the TLV. Suder et al.[21, 22] showed that the TLV interacts with the shock system, causes notable loss and blockage through mixing and diffusion. As in sub-figure (d), the TLV reduces the relative Mach number on its path before the bow shock. After the TLV-shock interaction, a tremendous amount of low momentum fluid is present. The above four span positions have been marked out with dashed horizontal lines in Fig. 6, the adjoint sensitivity plot. With the cost function defined as entropy generation, the adjoint gradient distribution on SS matches well with the position and strength of the shock system. There are also some other sensitive regions marked out by the four dashed lines, such as region A on PS, which corresponds to the passage shock on PS, see sub-figure (b) in Fig. 7. Besides the shock system, the tip leakage flow is also a big source for entropy generation. Historic literature[19, 20] claimed that the pressure difference between the SS and PS plays an important role in driving fluid through the tip gap. As demonstrated in Fig. 8, the TLV mainly initiates near LE and is continuously strengthened by jetlike fluids going over the blade tip until at least 70% chord length. Here the stream traces start from the middle grid line of the tip gap mesh blocks. The static pressure contours on the blade surface is also presented. On the PS near the blade tip, the flow becomes subsonic after the LE bow shock, as shown in the plot (d) of Fig. 7. Adjoint sensitivity says that this place needs to be thickened (see

5. Results and Discussions Here we present the adjoint sensitivities and design optimization results of the NASA Rotor 67 blade row. 5.1. Adjoint Sensitivity Analysis Figure 6 shows the adjoint sensitivity contours on the whole blade surface of the NASA Rotor 67 transonic fan blade. The blade surface is now cut open from the trailing edge (TE). The dashed vertical line in the middle of this figure is the leading edge (LE). The left half plot represents the pressure surface (PS) while the right half for the suction surface (SS). The x and y coordinates are grid nodes indexes on the blade surface. There are totally 145 × 65 RBFs, equal to the number of grid nodes on the blade surface. We place RBFs exactly on grid nodes such that a good resolution of adjoint sensitivity can be expected. In order to resolve the flow near the blade tip and casing, grid nodes cluster there. In the radial direction, the 25th , 40th , 44th , 50th and 58th RBF locate at 17.6%, 73.4%, 86.0%, 95.8% and 99.3% span, respectively. There is similar cluster effect in the chordwise direction. The 60th , 100th and 120th RBF correspond to 9.6% chord length on PS, 22.8% on SS and 65.4% on SS. All these locations have been marked by the thin dashed lines in Fig. 6. 7

4.0E-07 3.5E-07 3.0E-07 2.5E-07 2.0E-07 1.5E-07 1.0E-07 5.0E-08 0.0E+00 -5.0E-08 -1.0E-07 -1.5E-07 -2.0E-07 -2.5E-07 -3.0E-07 -3.5E-07 -4.0E-07

TE Tip 60 50

PS 20

60 C

22.8%

80

100

SS

65.4%

120

TE 140

B

60

99.3% 95.8%

50 A

86.0%

40

LE

9.6%

40

73.4%

40

30

30 17.6%

20

20

10

10

Hub

20

40

60

80

100

120

140

Figure 6: Adjoint sensitivity

Ma: 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 Y

Z

X

SS

0.9 PS

1.3

4 1.

1.0

SS

1.0 0.9 1.1

PS

(a) 73.4% Span

0.5

Figure 8: Tip clearance flow

1.4

1.3 1.2

1.0 1.1

(b) 86.0% Span

(c) 95.8% Span

Table 2 shows the CPU times in each step of one optimization cycle, normalized by the time of flow solution. The time costs of a brute force mesh deformation and the localized two-level mesh deformation are compared. Note that even in the “brute force” mesh deformation, some simplifications can still be adopted. For example, based on the fast decaying property of RBF bumps, we can ignore the influences of surface nodes that are far away enough when determining the movements of nodes on a certain mesh line. For the “brute force” case presented in Tab 2, the influences of surface nodes that are more than 5

(d) 99.3% Span

Figure 7: B2B relative Mach number contours near peak efficiency

region B and C marked by green ellipses in Fig. 6). The explanation is that thickening of B and C reduces the area of local flow passage such that the velocity of the fluid is increased, resulting in lower static pressure. This obviously decreases the pressure difference between PS and SS, thus suppresses the tip clearance flow.

Table 2: Normalized computation time

5.2. Optimization Procedure The optimization is conducted in a parallel environment using 6 CPUs (3.6GHz each). The total 430,000 mesh nodes are divided into 11 blocks. Heaviest computation load for a single CPU is 21%.

Brute Local

8

Flow 1.000 1.000

Adjoint 2.074 2.074

Gradient 2.486 0.222

New Mesh 0.043 0.003

Total 5.603 3.299

−4

sgen or I

1.9

x 10

1.85

0.938

1.8

0.936

1.75

0.934

Tip

Tip

0.940

sgen I η

TE

1.7 0

1

2

3

4

5 6 7 Design Cycle

8

9

10

11

LE

TE

LE

0.932 12 Hub

Hub

Figure 9: Adjoint optimization process

Ref Opt

−4

I 10 1.858 1.782





−4

sgen 10 1.858 1.732



m ˙ (kg/s) 33.52 33.76

(b) PS

Figure 10: Geometries comparison: blue-reference; green-optimized.

Table 3: Overall performances

η (%) 93.32 93.82

π 1.625 1.627

1.4 At 73.4% Span

Isentropic Mach



(a) SS

half decay radius away are ignored when calculating the movements of nodes on a certain mesh line. Also, the mesh deformation is limited in the O mesh surrounding the blade for both the “brute” and local cases. This is because the deformation magnitudes have already become negligible when they reach the outer boundary of the O mesh. The localized two-level method significantly cuts down the CPU time in mesh deformation. As demonstrated in Tab. 2, the adjoint solution takes about twice the time of flow solution because, on the one hand, the evaluation of terms in the adjoint equations takes longer time than that in flow equations, on the other hand, one-third more multi-grid iterations are used in the adjoint solution to ensure good convergence. The gradient calculation step takes much longer time than the new mesh generation not only for a two order higher operational count in mesh deformation (see Tab. 1) but also for the gradient collection procedures for every RBFs. Perturbations of 3 sub-blocks of the O mesh surrounding the blade are implemented parallelly. Then, perturbation of the tip gap block begins, with its boundary nodes movements equal to the movements of the O mesh layer surrounding the blade tip. As shown in Fig. 9 and Tab. 3, the cost function I and entropy generation sgen drop by 4.1% and 6.8% respectively after 10 design cycles. The adiabatic efficiency is optimized from 93.32% to 93.82% while the mass flow rate and total pressure ratio increase by 0.72% and 0.12%.

1.2

Reference Optimized

1.0

0.8

0.6 0.0

0.2

0.4

0.6

0.8

1.0

x/c (a) 73.4% Span 1.4

Isentropic Mach

At 99.3% Span

1.2

1.0

0.8

0.6 0.0

Reference Optimized

0.2

0.4

0.6

0.8

1.0

x/c (b) 99.3% Span Figure 11: Isentropic Mach number before/after optimization

lines represent the old geometry and green for the optimized. On SS (plot (a)), the geometry modifications near the shock can be seen clearly. For PS (plot (b)), the blade tip region is thickened. These geometric variations are directly decided by the adjoint gradient on the blade surface, see Fig. 6. The isentropic Mach number on the blade surfaces at different span positions before and after optimization are

5.3. Details of Design Optimization Figure 10 compares the blade geometry of Rotor 67 before and after design optimization. Here the blue mesh 9

Figure 13: Flow separation on SS before/after optimization

(a) Reference

(b) Optimized

blown downstream now maintains its circumferential momentum. The changes of stream traces can only be caused by the diminishing of the TLV, in both strength and range of influence. Besides the tip clearance flow, optimization of flow separation (caused by shock-boundary interaction) on SS is also investigated. Though three-dimensional flow separation is complicated, especially in turbomachinery, Gbadedo et al.[23] testified that it can be predicted by RANS codes with good accuracy, particularly the limiting streamline patterns. Fig. 13 provides limiting streamlines on SS together with contours of static pressure. As studied by Arnone[24], for NASA Rotor 67 blade near peak efficiency, the passage shock and separation line coincide on the outer span of SS. This flow separation is not severe and reattaches quickly. In our study, the separation line in the reference flow field begins near the blade tip and stretches down to 70% span. After optimization, the weaker shock causes a smaller separation, which ends at about 80% span. The spanwise distributions of total pressure ratio, total temperature ratio and flow angle at the exit are examined to better manifest the aerodynamic improvements. As demonstrated in Fig. 14(a), the optimized blade row increases the total pressure ratio in above 30% span. However, it decreases the total pressure ratio from 10% to 30% span. This total pressure defect is attributed to the enlargement of flow separation near the trailing edge, see Fig. 13. Fig. 14(b) shows that the work done to fluid per unit mass flow above 60% span is slightly less in the optimized case. This helps to achieve higher adiabatic efficiency. The last sub-figure gives the flow turning angle in degrees. In the outer half span, the flow turning angle becomes smaller than the original Rotor 67. The decrease in flow turning agrees with the reduced work done to the fluid per unit mass flow. The mechanical/structural design is equally important as the aerodynamic design. The spanwise distributions of

Figure 12: Tip clearance flow before/after optimization

shown in Fig. 11. It is clear that our optimization weakens and delays the shock on SS. Explanations have been given in the adjoint sensitivity analysis part. To demonstrate how the design optimization influences the tip clearance flow, we examine the stream traces starting at the middle line of the tip gap mesh block in Fig. 12, together with the static pressure contours. As marked out by the dash-dotted blue circles, in the optimized design, the pressure at the blade tip near the PS leading edge is notably lower than the reference value. This reduces the pressure difference between PS and SS where the tip clearance flow mainly initiates, thus a weaker tip leakage vortex (TLV) can be expected. After the leading edge, jet-like flow through the tip gap along the chord would also be weaker since the optimized geometry thickens the blade tip on PS, see Fig. 10(b). The two factors together suppress the tip clearance flow, resulting in a weaker TLV. Since both the shock system and the TLV are weakened, their interaction becomes less harmful. The corresponding flow diffusion and blockage would be relieved. The above analysis can be proved by comparing the stream traces in the tip region before and after optimization. In Fig. 12, all the stream traces start in the tip gap of blade 1 and go towards blade 2 and 3. After the design optimization, part of the fluid which originally rolls into the TLV now maintains its circumferential momentum after traveling across the passage between blade 1 and 2, then go over the tip gap of blade 2 into the next flow passage, as marked out by the dashed blue circle at the middle chord of blade 2. Similar actions in the passage between blade 2 and 3 are observed and marked out by the solid blue circle, where the fluid originally caught by the TLV and 10

1.0

1.0

0.8

0.8

0.8

0.6

0.6

0.6

r

r

r

1.0

0.4

0.4

0.2

0.2

0.0 1.60

1.62

π

1.64

(a) Total pressure ratio

1.66

0.0 1.11

0.4

Reference Optimized 1.16

θ

1.21

1.26

0.2 0.0 −55

(b) Total temperature ratio

−50

−45

β

−40

−35

−30

(c) Flow turning angle (deg)

Figure 14: Spanwise distribution of π, θ and β at outlet

the blade maximum thickness and its chordwise location before and after design optimization are shown in Fig. 15. The left plot shows the maximum thickness location. The middle plot shows the maximum thickness while the right plot is a zoom-in of region A in the middle plot. Generally, the maximum thickness of the blade profile does not decrease in the current design optimization. Though the location of the maximum thickness does change. It moves upstream in 50% to 70% span and downstream in 70% to 100% span. Usually, in a rotor blade without sweep and lean, due to the centrifugal force, the maximum stress happens near the stacking line below 50% span. Careful Finite Element Analysis (FEA) steps are required to ensure structural/mechanical feasibility of a new blade shape. The focus of the present paper is on the aerodynamics. However, based on observations, some discussions are listed here. 1. In the current design, the maximum thickness and its location are not changed below 50% span. 2. The pressure surface redesign near blade tip does not change the size of the tip clearance, i.e., the height of the tip gap. 3. As indicated by the maximum thickness change in plots (b) and (c) in Fig. 15, the geometry change on PS near blade tip is small. The maximum thickness at the tip is increased by only 1/30 of the original value. And the change is limited in 2 percent of the blade height. Thus, the modification of blade weight due to PS redesign near blade tip is small. 4. The PS redesign near blade tip can not make the maximum stress worse since the stress near blade tip is small compared with the largest stress on the blade. 5. The PS redesign near blade tip might not vary the vibration frequencies much, since the corresponding change in blade weight is small. Considering the discussions above, the authors believe that the geometric change on the PS near the tip is feasible in the structural sense.

11

6. Conclusion An optimization platform for turbomachinery blades with complex mesh configuration in a parallel computation environment is built. Based on the same parallel framework of the flow solver, we code a continuous adjoint solver for three-dimensional viscous flows in rotors. To ensure convergence of the adjoint field, a non-reflection inlet/outlet boundary condition is adopted. A localized two-level mesh deformation method is designed to cut down the computation cost in grid perturbation while providing reliable meshes for both adjoint gradient collection and next cycle CFD calculation. With a Gaussian radial basis function (RBF) as basic bump, this method works well not only for the O mesh surrounding the blade but also for the O-H mesh blocks inside the tip gap. The adjoint sensitivities on the whole blade surface of the NASA Rotor 67 are presented and discussed in detail. Bridges between the adjoint sensitivities and physical phenomena are constructed. Aerodynamic optimization for high adiabatic efficiency is conducted for this blade row with a mass flow rate constraint. After the redesign of the blade profile, the strength of the shock system is reduced. The tip leakage vortex (TLV) is notably suppressed through small changes of the pressure surface geometry near the blade tip. The weakening of shock and the TLV together relieve the flow diffusion and blockage related to shock-TLV interaction. Another advantage of shock strength reduction is that the flow separation caused by the passage shock-boundary layer interaction on the suction surface becomes weaker. Future work contains four parts. Firstly, practical constraints like maximum thickness, total pressure ratio, and blade weight will be added. Secondly, multi-disciplinary optimizations will be explored, including structural stress and vibration. The third target is multipoint and multiobjective optimizations. Fourthly, detailed investigations into the adjoint field could be conducted (around shock wave and tip clearance) so that the physics behind it might be better understood.

1.0

1.0

1.00

A 0.8

0.98

0.6

0.6

0.96

0.4

0.4 0.2

r

r

r

0.8

0.2 Reference Optimized

0.0 40.0 50.0 60.0 70.0 80.0 Maximum thickness location (percent of local chord)

(a) Maximum Thickness Location

0.94 0.92

Reference Optimized

0.0 0.0 2.0 4.0 6.0 8.0 10.0 Maximum thickness (percent of local chord)

Reference Optimized 0.90 2.0 2.5 3.0 3.5 4.0 Maximum thickness (percent of local chord)

(b) Maximum Thickness

(c) Region A

Figure 15: Maximum thickness and its location before/after optimization

Acknowledgements

[12] Yang, S., Wu, H., and Liu, F., “Aerodynamic Design of Cascades by Using an Adjoint Equation Method,” AIAA Paper 2003-1068, 2003. [13] Wu, H. and Liu, F., “Aerodynamic Design of Turbine Blades Using an Adjoint Equation Method,” AIAA Paper 2005-1006, 2005. [14] Walther, B. and Nadarajah, S., “Adjoint-Based Constrained Aerodynamic Shape Optimization for Multistage Turbomachines,” Journal of Propulsion and Power , Vol. 31, No. 5, September-October 2015. [15] Strazisar, A. J., Wood, J. R., Hathaway, M. D., and Suder, K. I., “Laser Anemometer Measurements in a Transonic Axial-Flow Fan Rotor,” Tech. Rep. 2879, NASA, 1989. [16] Jameson, A., Schmidt, W., and Turkel, E., “Numerical Solution of the Euler Equations by Finite Volume Methods Using RungeKutta Time-Stepping Schemes,” AIAA Paper 81-1259, 1981. [17] Menter, F. R., “Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications,” AIAA Journal, Vol. 32, No. 8, 1994, pp. 1598–1605. [18] Strazisar, A. J., “Investigation of Flow Phenomena in a Transonic Fan Rotor Using Laser Anemometry,” ASME Journal of Engineering for Gas Turbine and Power , Vol. 107, No. 2, April 1985, pp. 427–435. [19] Pandya, A. and Lakshminarayana, B., “Investigation of the Tip Clearance Flow Inside and at the Exit of a Compressor Rotor Passage. 1. Mean Velocity-Field,” ASME Journal of Engineering for Power , Vol. 105, No. 1, 1983, pp. 1–12. [20] Dawes, W. N., “A Numerical Analysis of the Three-Dimensional Viscous Flow in a Transonic Compressor Rotor and Comparison With Experiment,” ASME Journal of Turbomachinery, Vol. 109, No. 1, Jan 1987, pp. 83–90. [21] Suder, K. L. and Celestina, M. L., “Experimental and Computational Investigation of the Tip Clearance Flow in a Transonic Axial Compressor Rotor,” ASME Journal of Turbomachinery, Vol. 118, No. 2, April 1996, pp. 218–229. [22] Suder, K. L., “Blockage Development in a Transonic, Axial Compressor Rotor,” ASME Journal of Turbomachinery, Vol. 120, No. 3, July 1998, pp. 465–476. [23] Gbadebo, S. A., Cumpsty, N. A., and Hynes, T. P., “ThreeDimensional Separations in Axial Compressors,” ASME Journal of Turbomachinery, Vol. 127, No. 2, April 2005, pp. 331– 339. [24] Arnone, A., “Viscous Analysis of 3-Dimensional Rotor Flow Using a Multigrid Method,” ASME Journal of Turbomachinery, Vol. 116, No. 3, July 1994, pp. 435–445.

The authors would like to thank The National Natural Science Foundation of China (Grant No. 51376009 and No. 51676003) for supporting the present research work. References [1] Hicks, R. and Henne, P., “Wing Design by Numerical Optimization,” Journal of Aircraft, Vol. 15, No. 7, 1978, pp. 407–413. [2] Luo, J., Xiong, J., Liu, F., and McBean, I., “Three-Dimensional Aerodynamic Design Optimization of a Turbine Blade by Using an Adjoint Method,” ASME Journal of Turbomachinery, Vol. 133, No. 1, 2011. [3] Luo, J., Zhou, C., and Liu, F., “Multi-Point Design Optimization of a Transonic Compressor Blade by Using an Adjoint Method,” ASME Journal of Turbomachinery, Vol. 136, No. 5, 2014. [4] Luo, J., Liu, F., and McBean, I., “Turbine Blade Row Optimization Through Endwall Contouring by an Adjoint Method,” Journal of Propulsion and Power , Vol. 31, No. 2, March-April 2015. [5] Marta, A. C., Shankaran, S., Wang, Q. Q., and Venugopal, P., “Interpretation of Adjoint Solutions for Turbomachinery Flows,” AIAA Journal, Vol. 51, No. 7, Jul 2013, pp. 1733–1744. [6] Jameson, A., “Optimum Aerodynamic Design Using CFD and Control Theory,” AIAA Paper 1995-1729, 1995. [7] Lee, B. J. and Kim, C., “Automated Design Method of Turbulent Internal Flow Using Discrete Adjoint Formulation,” Aerospace Science and Technology, Vol. 11, No. 2-3, MAR-APR 2007, pp. 163–173. [8] Walther, B. and Nadarajah, S., “Optimum Shape Design for Multirow Turbomachinery Configurations Using a Discrete Adjoint Approach and an Efficient Radial Basis Function Deformation Scheme for Complex Multiblock Grids,” ASME Journal of Turbomachinery, Vol. 137, No. 8, Aug 2015. [9] He, L. and Wang, D. X., “Concurrent Blade Aerodynamic-Aeroelastic Design Optimization Using Adjoint Method,” ASME Journal of Turbomachinery, Vol. 133, No. 1, Jan 2011. [10] Huang, H. and Ekici, K., “A Discrete Adjoint Harmonic Balance Method for Turbomachinery Shape Optimization,” Aerospace Science and Technology, Vol. 39, DEC 2014, pp. 481–490. [11] Jameson, A., “Aerodynamic Design via Control Theory,” Journal of Scientific Computing, Vol. 3, No. 3, June 1988.

12