Aerodynamic shape optimization of hovering rotor blades using a Non-Linear Frequency Domain approach

Computers & Fluids 51 (2011) 1–15

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Aerodynamic shape optimization of hovering rotor blades using a Non-Linear Frequency Domain approach Charles A. Tatossian 1, Siva K. Nadarajah ⇑,2, Patrice Castonguay Department of Mechanical Engineering, McGill University, Montreal, Quebec, Canada

a r t i c l e

i n f o

Article history: Received 9 December 2009 Received in revised form 13 September 2010 Accepted 23 June 2011 Available online 27 July 2011 Keywords: Adjoint approach Optimization Non-Linear Frequency Domain method Rotorcraft

a b s t r a c t This paper presents a discrete adjoint-based aerodynamic optimization algorithm for helicopter rotor blades in hover using a Non-Linear Frequency Domain approach. The procedure to be used can be broken down into two phases: first, investigate and develop an accurate inviscid numerical model for the current generation of helicopter rotor blades; and second, demonstrate the approach for the redesign of a transonic Caradonna–Tung two-bladed rotor. The results in determining the optimum aerodynamic configurations require an objective function which minimizes the inviscid torque coefficient and maintains the desired thrust level at transonic conditions. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Rotorcraft aerodynamicists continually challenge themselves to improve the aerodynamic performance of today’s helicopter blade designs. Modern designs in applied aerodynamics can be resolved via two fundamental approaches: theoretical fluid dynamics or computational fluid dynamics. Although both methodologies have been vigorously scrutinized to study the interaction between the fluid flow and solid bodies, computational fluid dynamics has always been characterized as a supplemental aid in the design process of rotorcraft rather than serving as a direct tool. Despite the recent efforts in ameliorating existing steady flow aerodynamic shape optimization algorithms, there remains a considerable need to develop innovative and cost-effective optimal control techniques for the design of aerodynamic surfaces subjected to unsteady loads. Researchers have applied aerodynamic shape optimization algorithms to numerous steady-state problems, ranging from the design of two-dimensional airfoils to full aircraft configurations to decrease drag, increase range, promote lift, etc. [1–11]. However, unlike fixed wing aircraft, helicopter rotors operate in an unsteady flow and are constantly subjected to unsteady loads. Reliable

⇑ Corresponding author. 1 2

E-mail address: [email protected] (S.K. Nadarajah). Graduate Student. Associate Professor Student.

0045-7930/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2011.06.014

predictions of helicopter rotors are challenging, primarily due to the heavy dependency on the resolution of the blade-vortex interaction near the tip region of the blade. This interaction strongly influences the incoming flow and consequently alters the effective angle of attack. As the helicopter moves forward, the rotor blade comes across a combination of rotational and translational flow, and the symmetry formed during hover vanishes. Local supersonic zones that terminate at a shock wave develop on the advancing side, while at the retreating phase, the blade’s velocity relative to the air decreases and the blade approaches the stall angle. This causes significant flow separation to occur on the upper surface of the blade which in turn produces a loss in lift. All these issues have to be carefully taken into consideration to fully comprehend the performance of helicopter rotor blades and to truly appreciate the simplicity of the final design. Despite these challenges, numerical simulation of the flow over helicopter rotors have been attempted and perfected by a number of researchers such as Srinivasan et al. [12], Ahmad and Strawn [13], Datta et al. [14], Hariharan and Sanker [15], Yeshala et al. [16], and Pomin and Wagner [17]. Leishman [18] presents an extensive list of papers, reports, and articles related to computational methods for helicopter aerodynamics. Presently, a very limited number of research work for aerodynamic shape optimization of helicopter blades exists. The transition of adjoint-based algorithms originally developed for fixed-wing towards three-dimensional unsteady viscous flows has been slow primarily due to the demanding computational cost. Nevertheless,

2

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15

Nomenclature D E F G I M R b R S S V b f p

artificial dissipation flux internal energy numerical flux vector gradient cost function Mach number residual Fourier coefficient of residual shape function face areas of computational cell cell volume boundary velocity component flux vector pressure

Nadarajah and Jameson [19,20] have pursued the development of optimum shape design for two- and three-dimensional unsteady flows using the time accurate adjoint based design approach. Nadarajah derived and applied the time accurate adjoint equations (both the continuous and discrete) to the redesign of an oscillating airfoil in an inviscid transonic flow. The redesigned shape achieved a reduction in the time-averaged drag while maintaining the timeaveraged lift. The approach utilized a dual time stepping [21] technique that implements a fully implicit second-order backward difference formula to discretize the time derivative. Encouraging results were obtained at a substantial computational expense. In the case of helicopter rotors, Lee et al. [22], presented the redesign of the Caradonna–Tung and UH-60 rotors via a continuous adjoint approach, along with a solution-adaptive mesh refinement method to efficiently capture the blade’s tip vortex while simultaneously reducing computational cost. In 2007, Morris et al. [23] provided a generic optimization tool which allowed high-fidelity aerodynamic optimization of two- and three-dimensional blades. Heavily constrained two-dimensional cross-sections showed significant improvements over existing designs. Despite these efforts, the optimization of unsteady flows continues to present severe challenges in terms of computational cost. As a result, there has been much effort focused on the development of efficient and practical alternatives to the study of unsteady periodic problems. One potential approach is to use periodic methods. Linearized frequency domain and deterministic stress methods [24] are examples of periodic methods that are widely used in industry. Unfortunately, the inability of these methods to accurately model the solution becomes evident for systems that contain strong nonlinearities. The Harmonic Balance technique, a pseudospectral approach initially proposed by Hall [25] and later modified by McMullen [26,27], has been validated against both the Euler and Navier–Stokes equations for a number of unsteady periodic problems and has been shown to account for strong nonlinearities. The cost associated with spectral methods like McMullen’s NonLinear Frequency Domain method (NLFD) is proportional to the cost of the steady state solution multiplied by the number of desired temporal modes. For inviscid flow, McMullen [28] has shown that to accurately model an oscillating airfoil pitching about its quarter chord, a temporal resolution of only 1 mode above the fundamental frequency (or equivalently three time samples per period) is needed using the NLFD method versus the 45 time samples needed with a backward difference formulation of the time derivative [21]. These results demonstrate the potential of the method to provide significant reduction in computational cost for the analysis and design of more realistic problems such as helicopter rotors, turbomachinery, and other unsteady devices operating

t tf u w ^ w x

a n

q w xr i, j, k k

time final time velocity (physical domain) state vector Fourier coefficient of state vector coordinates (physical domain) angle of attack coordinates (computational domain) density Lagrange multiplier reduced frequency cell indices wave number

in the transonic regime. Nadarajah et al. [29] extended their twodimensional optimum shape design for unsteady flows from a time accurate scheme to the NLFD approach and in the process developed the NLFD adjoint equations. The method was further extended for three-dimensional inviscid and viscous flows and a complete comparison between the two techniques [30,31] was presented. In order to quantify the required number of time steps per period required for both the time accurate and the NLFD approaches, a temporal resolution study was performed in [32], where the time accurate code was run with 5, 11, 23, 47, and 767 time steps per period and 3, 5, and 7 for the NLFD approach for a model two-dimensional inviscid test case. The 767 time step case was used as the control solution and the errors in the time averaged lift coefficient for all test cases for both the time accurate and NLFD methods were compared to that given by the control solution. A RAE airfoil on a 193  33 grid was employed and the computations were performed at a freestream Mach number of 0.78, at a mean angle of attack, ao = 0o, a maximum angle of attack, am = 1.01o, and at a reduced frequency, xr = 0.20. Nadarajah [32] demonstrated that the error in the time-averaged life coefficient as well as the accuracy of the gradient computed via the adjoint approach exhibited a spectral-type convergence for the NLFD method when compared to second-order accuracy for the timeaccurate method. The numerical simulation of flow over helicopter rotor blades via a frequency domain based method is currently being pursued by a number of researches. Choi et al. [33] presented a time-spectral approach to simulate the UH-60A rotor. The conclusion drawn from the research suggested that the time-spectral method provides great potential for adjoint-based optimization algorithms for periodic problems when compared to time-accurate schemes. Kumar and Murthy [34] introduced a multiblade coordinate system (MCS) that offers the potential to further reduce the minimum number of harmonics around the azimuth to accurately capture the flow around the helicopter blade. Finally, Butsuntorn and Jameson [35] recently proposed a time-spectral method together with a vorticity confinement formulation to simulate the unsteady threedimensional rotorcraft flow for both the Euler and Reynolds averaged Navier–Stokes equations. The modeling of unsteady aerodynamic design sensitivities using either the harmonic balance technique or the Non-Linear Frequency Domain approach have also been investigated by Duta et al. [36] and Thomas et al. [37]. Duta et al. [36] have presented a harmonic adjoint approach for unsteady turbomachinery design. The aim of the work was to reduce blade vibrations due to flow unsteadiness. The research produced adjoint methods that were based on a linearized analysis of periodic unsteady flows. Thomas

3

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15

et al. [37] presented a viscous discrete adjoint approach for computing unsteady aerodynamic design sensitivities. The adjoint code was generated from the harmonic balance flow solver with the use of an automatic differentiation software compiler. The goal of this work is to extend our three-dimensional NLFD adjoint-based design framework for helicopter rotor blades in hover, and investigate its impact by analyzing the overall change in aerodynamic performance. To simulate and redesign the flow over helicopter rotors, a number of additional developments have been added to the present framework such as hover boundary conditions, periodic boundary conditions and new adjoint boundary conditions. The flow solver is initially validated for the Caradonna–Tung experimental rotor blade at a tip Mach number of 0.520 and collective pitch of 0° and lifting cases at a tip Mach number of 0.439 and 0.877 and collective pitch of 8° and the blade is subsequently redesigned using the NLFD adjoint-based algorithm. The number of modes required for the accurate representation of the flow solution and the gradient are explored. The motivation of the research has been fueled both by the success of our current capability for automatic shape optimization for unsteady flows and the future potential of the NLFD method.

represents the projection of the ni cell face along the xj axis. These equations can be written in integral form for a domain D with boundary B as

2. Governing equations and their discretization

2.1. Spatial discretization

The Cartesian coordinates and velocity components are denoted by x1, x2, x3, and u1, u2, and u3. Einstein notation simplifies the presentation of the equations, where summation over i = 1–3 is implied by a repeated index i. The three-dimensional Euler equations then take the form,

The convective flux is represented in discrete form for each computational cell using a central second-order discretization. Eq. (1) can then be written for each computational cell as

@w @fi ¼ 0 in D; þ @t @xi

The residual can then be represented as

where the state vector w and inviscid flux vector f are described respectively by

9 8 q > > > > > > > > > > > = < qu1 > w ¼ qu2 ; > > > > > > > qu3 > > > > > ; : qE

9 8 qðui  bi Þ > > > > > > > > > > > = < qu1 ðui  bi Þ þ pdi1 > and f i ¼ qu2 ðui  bi Þ þ pdi2 : > > > > > > > qu3 ðui  bi Þ þ pdi3 > > > > > ; : qEðui  bi Þ þ pui

In these definitions, q is the density, bi are the Cartesian velocity components of the boundary, E is the total energy and dij is the Kronecker delta function. The pressure is determined by the equation of state



 1 p ¼ ðc  1Þq E  ðui ui Þ ; 2 and the stagnation enthalpy is given by

H ¼Eþ

p

q



Z

w dV þ

Z

D

fj dAi ¼ 0;

B

where dAi is the component of area projected in the xi direction. A finite volume scheme is derived by applying the equation above directly to control volumes to give a set of ordinary differential equations of the form

dðwVÞ þ RðwÞ ¼ 0 in D; dt

 @xi ; @nj

V i;j;k

@wi;j;k þ RðwÞi;j;k ¼ 0: @t

K 1 ij ¼



 @ni : @xj

The Euler equations can then be written in computational space as

@ðJwÞ @F i ¼ 0 in D; þ @t @ni where the inviscid flux contribution can be defined with respect to the computational cell faces by Fi = Sijfj, and the quantity Sij ¼ JK 1 ij

ð2Þ

RðwÞi;j;k ¼ hiþ1;j;k  hi1;j;k þ hi;jþ1;k  hi;j1;k þ hi;j;kþ1  hi;j;k1 ; 2

2

2

2

2

2

ð3Þ

where hiþ1;j;k ¼ F iþ1;j;k  Diþ1;j;k and D represents the artificial dissipa2 2 2 tion term. The  12 notation indicates that the quantity is calculated at the flux faces. The values of the flow variables are stored at the cell centers and can be regarded as cell averages. Accordingly, the convective flux F iþ1;j;k at the cell face is computed by taking the aver2 age of the flux contributions from each cell across the cell face. A CUSP scheme is employed for the artificial dissipation [38]. The spectral radius of the flux Jacobian matrix are rescaled in the each direction to better handle the high aspect ratio cells in the vicinity of the blade tip. 2.2. Temporal discretization The derivation of the NLFD method starts with Eq. (1) and assumes that the vector of flow variables w and local residual R can be represented by separate Fourier series:

w¼

J ¼ detðKÞ;

ð1Þ

where V is the cell volume, and the residual R(w) is evaluated by summing the fluxes through the cell faces. The simulations contained in this research are restricted to rigid mesh translation; therefore, the cell volumes remain constant in time. To eliminate odd–even decoupling of the solution and overshoots before and after shock waves, the convective flux is added to a diffusion flux. The artificial dissipation scheme used in this research is based on a convective upwind and split pressure (CUSP) scheme.

N

;

where c is the ratio of the specific heats. For discussion of real applications using a discretization on a body conforming structured mesh, it is also useful to consider a transformation to the computational coordinates (n1, n2, n3) defined by the metrics

K ij ¼

@ @t

1 2 X

N

^ k eikt ; w

R¼

k¼N2

1 2 X

b k eikt R

ð4Þ

k¼N2

pffiffiffiffiffiffiffi where i ¼ 1. The Fourier representations are then substituted into the semi-discrete form of the governing equations as described in Eq. (1) to yield,

2N 3 N 1 1 2 2 X d 4X ikt 5 b k eikt ¼ 0 in D: ^ þ R V wk e dt N N k¼ 2

ð5Þ

k¼ 2

As a result,

c bk; ^k þ R R k ¼ ikV w

ð6Þ

forms the new unsteady residual in the frequency domain for each wavenumber and must be solved iteratively. The solver attempts to

4

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15

find a solution, w, that drives this system of equations to zero for all wavenumbers. The nonlinearity of the unsteady residual and, hence, the bulk of our computing effort, stems from the spatial operator, R. There are two approaches to calculating the spatial operator expressed in the frequency domain. The first uses a comb k directly from w ^ k. plex series of convolution sums to calculate R This approach is discarded due to its massive complexity (considering artificial dissipation schemes and turbulence modeling) and cost that scales quadratically with the number of modes N. This paper implements a pseudo-spectral approach that relies on the computational efficiency of the Fast Fourier Transform (FFT). To calculate R in the frequency domain, several transformations between the physical and frequency domains are performed by FFT. The computational cost of this transform scales like N log (N), a significant savings over a similar method that uses convolution sums and scales like N2. A diagram detailing the transformations used by the pseudo spectral approach at each stage of the modified multistage Runge– Kutta scheme is provided in Fig. 1. The pseudo-spectral approach begins by initializing the state vector, w(t), at all time instances. At the first iteration, w(t) will take on the value of the initial condition and for subsequent iterations, the values are based on the previous iterations. At each of these time instances the steady-state operator R(w(t)) can be computed by summing the convective and artificial dissipation fluxes as mentioned in the previous subsection. A FFT is then used to transform the state vector and spatial operator to the frequency b k are known for all wavenumbers. The un^ k and R domain, where w b k to the specsteady residual c R k can then be calculated by adding R ^ k . The iteration tral representation of the temporal derivative ikV w ^ k is then ^ k is updated. Using an inverse FFT, w is advanced and the w transformed back to the physical space resulting in a state vector w(t) sampled at evenly distributed intervals over the time period. The wall and far-field boundary conditions are imposed within the time domain. Consistent with the time accurate approach, we can numerically integrate our residual in fictitious time t⁄ resulting in the following equation:

V

the memory cost increases. Second, ’calls’ to FFT routines must be implemented within the time-stepping scheme to perform both the FFT and inverse FFT. An important limitation of the current implementation is the assumption that the cell volume, V, is constant in time. In the work presented in this paper, a rigid body translation and rotation is employed, and, therefore, this assumption is valid. An extension of the method for unsteady motions, where small grid deformations are present, the volume, V(t), will be a function of time and must be included in Eq. (5) within the time derivative. However, for cases, where large grid deformations are observed that might include grid adaptation which potentially could increase or decrease the number of cells to improve the quality of the grid for the flow solver, an interpolation of the state vector, w(t), to a common fixed point in space might be needed to transform the state vector to the frequency domain. Nevertheless, such an approach would reduce the order of the method and cast the approach as a reduced order model. The cases attempted in this work are limited to the simulation of rotor blades in a hover condition, where the collective pitch is fixed at all phases of the rotation. Therefore, an alternate cost effective approach to compute the flow would be to solve for the absolute velocities in the rotational frame of reference without physically rotating the computational grid such as that proposed by Holmes and Tong [39]. This approach requires the addition of a source term to account for the rotor angular velocity and is computational very attractive. However, it is our intention to further develop the technique and demonstrate the framework for rotor blades in forward flight, where the addition of a freestream Mach number as well as a variable collective pitch would require a full unsteady flow solver. The presented work is primarily to validate the approach for the hover flight condition.

^ k c dw  þ R k ¼ 0: dt

An unsteady residual exists for each wavenumber used in the solution and the pseudo-time derivative acts as a gradient to drive the absolute value of all of these components to zero simultaneously. A modified five stage Runge–Kutta time integration scheme is applied to march the solution to steady state. Local time stepping, implicit residual averaging, and multigrid are employed to accelerate the convergence. In synopsis, the NLFD approach can be implemented into any existing time-accurate flow solver with the following two primary modifications. First, the state vector must be modified to allow an additional dimension to hold the value at various time steps. This is unlike the time accurate approach, where the second-order backward difference approach only required the state vector to be stored at three time levels. For the NLFD approach if three or greater time steps are employed,

Fig. 1. Simplified data-flow diagram of the time advancement scheme illustrating the pseudo spectral approach used in calculating the non-linear spatial operator R.

Fig. 2. Periodic grid boundary.

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15

5

Fig. 3. Structured 257  65  49 C–H grid based on cylindrical topology for the Caradonna–Tung blade.

2.3. Grids and boundary conditions The computational grid employed for the validation study is a structured cylindrical C–H grid as illustrated in Figs. 2 and 3 for the Caradonna–Tung experiment [40]. The blade has a rectangular planform with an aspect ratio of six with NACA 0012 blade crosssections. The blade is untwisted with no precone angle. The MPI [41] domain topology is based on a N pi ¼ 4; N pj ¼ 1; N pk ¼ 3, where Np is the number of processors in each direction. Each processor contains a grid of size nx  nj  nk = 65  65  17. The total grid is 257  65  49. A concentration of grid points are placed at the blade’s leading and trailing edges as well as the blade tip to capture the tip vortex. Approximately half of the grids in the blade radial direction is placed within 20% of the blade tip. A near-periodic mesh is placed between blades as a one-to-one match at the periodic boundary was a challenge to generate with a C type mesh as illustrated in Fig. 2. Moreover, it would have forced undesirable high aspect ratio cells at the incoming periodic boundary, which would have both caused low grid quality as well as affected the convergence rate of the employed explicit-based solver. Fig. 2b demonstrates the concentration of points in the region of the leading edge forward of the blade tip. The inner radial mesh starts at the blade root; therefore, the blade root vortex was not modeled in this work. Every effort was made, however, to concentrate points at the blade tip within the limits of the mesh size. The outer radial shell as well as the top and bottom surfaces were placed at a distance of five rotor radii.

On the rotor surface, a flow tangency boundary condition is employed with a zeroth-order extrapolation of pressure across the wall. Periodic boundary conditions are imposed between blades in the blade azimuthal direction. Since non-congruent cells are used at the periodic boundary, an interpolation scheme based on the distance of the cell centers is used. Initially a preprocessing step is performed for the periodic boundary, where the cell center locations are computed and exchanged between processors. Then at the end of each multigrid cycle and only on the fine mesh, a periodic exchange of the state vector is executed. The state vector is then computed based on the ratio of the sum of product of the state vector in each cell within a certain radius and the inverse of its distance from the current cell and the sum of the inverse of the distance of all the cell centers from the current cell. The vector is then pre-multiplied by a transformation matrix to account for the rotation about the vertical axis. At the far-field boundary, which is composed of the upper, lower and outer-shell, hover boundary conditions have been implemented. The boundary conditions are based on the assumption that the implied presence of the rotor is to act as a three-dimensional point sink. The rotor outflow, which is situated underneath the blade, is represented via a one-dimensional momentum p theory ffiffiffi pffiffiffiffiffiffiffiffiffiffi equation, V e ¼ 2Mtip C t =2, and restricted to a radius of R= 2. On the other hand, the rotor inflow is characterized by a p point ffiffiffiffiffiffiffiffiffiffisink of mass, with its velocity modulus set to V r ¼ ðMtip =4Þ C t =2ðR=rÞ2 . The velocity, Vr, is a function of the tip Mach number, thrust coefficient, blade rotor radius, and radial distance from the center of the rotor disk. The velocity points radially inwards at each cell center at the outer boundaries and its magnitude is inversely proportional to the distance from the center of the rotor disk. These velocities were then used in the existing non-reflecting boundary conditions based on one-dimensional Riemann invariants normal to the boundary. The thrust coefficient, Ct, was set to the experimental value for each test case, as reported by Ahmad et al. [13], and its magnitude did not affect the final solution. However, it certainly improved the convergence rate of the flow solver as it prescribes a non-zero mass flux through the outer domain instead of an undisturbed flow at the far-field.

3. Optimization problem

Fig. 4. Design process.

The goal of this work is to minimize the torque generated by the rotor blade, by varying the shape at each blade section, subject to the constraint that the thrust coefficient is maintained. The objective is realized within this framework with the introduction of the adjoint approach for gradient-based optimization. There are generally two approaches to develop the adjoint equations: continuous or discrete. In the continuous adjoint approach the control theory is applied to the differential equations governing the flow. In the

6

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15

discrete adjoint approach, control theory is applied directly to the set of discrete field equations. A comparison of the two approaches were explored by Nadarajah et al. [42,43]. A discrete adjoint approach was employed for the work presented in this paper as described below as it provides for a more manageable way to derive the NLFD adjoint equations.





1 S1li1;j;k ATli;j;k ðwi;j;k  wi1;j;k Þ þ S1liþ1;j;k ATli;j;k ðwiþ1;j;k  wi;j;k Þ 2 2 2



T þ S2li;jþ1;k Ali;j;k ðwi;jþ1;k  wi;j;k Þ þ S2li;j1;k ATli;j;k ðwi;j;k  wi;j1;k Þ 2 2 



T þ S3li;j;kþ1 Ali;j;k ðwi;j;kþ1  wi;j;k Þ þ S3li;j;k1 ATli;j;k ðwi;j;k  wi;j;k1 Þ

RðwÞni;j;k ¼

2

3.1. Non-Linear Frequency Domain adjoint equations

2

The discrete adjoint equation is obtained by applying control theory directly to the set of discrete field equations. To formulate the discrete adjoint equation, we first define the cost function I as such

I¼

 N   2 1X -1 cnl  cnlinitial þ -2 cnd Dt; T n¼1

where cl and cd are the lift and drag coefficients at each time step, and N is defined as T = NDt, where Dt is the time discretization step. The next procedure in the development of the discrete adjoint equations is to differentiate the unsteady residual term, d c R k , from Eq. (6). There are two separate approaches to develop the NLFD discrete adjoint equations. In the first approach, we take a variation of the unsteady residual c R k represented in Eq. (6) with respect to the ^ state vector wk and shape function S, to produce

bk: ^ k þ dR dc R k ¼ ikVdw b k as a function of w ^ k . As The next step would be to expand d R mentioned earlier, this approach would require a series of convolub k as a function of dw ^ k . This method was not tion sums to express d R implemented due to its computational cost and added complexity. Instead, the adjoint equations were solved using a pseudo-spectral approach similar to the one applied to the flow equations. Here, the unsteady residual at each time step, RðwÞni;j;k , is differentiated to produce the discrete adjoint fluxes at each time step. The procedure is as follows; first, we take a variation of the residual term, which can be written as n

n

n

n

n

n

dRðwÞni;j;k ¼ dhiþ1;j;k  dhi1;j;k þ dhi;jþ1;k  dhi;j1;k þ dhi;j;kþ1  dhi;j;k1 ; 2

2

2

2

2

2

ð7Þ n dhiþ1;j;k 2

dF niþ1;j;k 2

dDniþ1;j;k . 2

with ¼  We then pre-multiply the variation of the discrete residual by the transpose of the Lagrange Multiplier and sum the product over the computational domain to produce the following jmax X imax X kmax N X 1X ðwn ÞT dRðwÞni;j;k Dt; T n¼1 i¼2 j¼2 k¼2 i;j;k

where imax, jmax, and kmax are the maximum dimensions in each coordinate direction, and dRðwÞni;j;k defines the variation of the residual of the time-dependent flow solution for n = 1, . . . , N. The adjoint equation can then be subsequently derived by the introduction of a Lagrangian function, where the variation of the cost function, dI, is augmented by the product of the variation of the discrete governing equation and the Lagrange Multiplier ðwni;j;k ÞT ,

dL ¼ dI 

jmax X imax X kmax N X 1X ðwn ÞT dRðwÞni;j;k Dt: T n¼1 i¼2 j¼2 k¼2 i;j;k

2

þDiþ1;j;k  Di1;j;k þ Di;jþ1;k  Di;j1;k þ Di;j;kþ1  Di;j;k1 ;

ð8Þ

2

2

2

2

2

ð9Þ

@f where AT is the transpose of the convective flux Jacobian @w , and Diþ1;j;k is the partial discrete adjoint artificial dissipation term. 2 To acquire the discrete adjoint boundary condition, Eq. (9) can be expanded along the cells, (i, 2, k), over the surface of the wing. In these cells, the terms from the objective function as well as the variation of the residual terms (7) that are multiplied to dwni;2;k are combined to form the discrete adjoint boundary condition. The boundary condition thus appears as a source term added to the total residual as shown below,

U ¼ S2;l wnli;2;k þ 

2-2 M21 p1 c

c

2-1 ½ðS2;2 sin a  S2;1 cos aÞ  clinitial 2 cM21 p1 c !

ðS2;1 sin a þ S2;2 cos aÞ

@p dwn : @w i;2;k

The values of the Lagrange multipliers for the cells below the wing are computed based on a zeroth-order extrapolation. Once the discrete adjoint fluxes and boundary conditions are established, the unsteady adjoint equations can be formulated. The semi-discrete from of the NLFD discrete adjoint equations can be expressed as

V

@w þ RðwÞ ¼ 0; @t

where R(w) is as defined in Eq. (9). Next, we assume that the adjoint variable and spatial operator can be expressed as a Fourier series: N

w¼

1 2 X

N

^ k eikt ; w

k¼N2

RðwÞ ¼

1 2 X

d k eikt : RðwÞ

ð10Þ

k¼N2

The derivation of the NLFD adjoint then follows that of the NLFD flow equations. The NLFD adjoint equations are expressed as

V

^k @w d  ¼ 0; þ RðwÞ k @s

d  ¼ ikV w d k . The pseudo-spectral approach illus^ k þ RðwÞ where RðwÞ k trated in Fig. 1 is employed in the NLFD adjoint code to form the unsteady residual. This term, in conjunction with a pseudo time derivative, provides an iterative solution process consistent with that documented for the flow equations. The gradient for the discrete adjoint is obtained by perturbing each mesh point on the rotor blade. Once the gradient G has been determined, it can be used to drive a variety of gradient-based search procedures. The search procedure used in this work is a descent method in which small steps are taken in the negative gradient direction. To accelerate the convergence, the gradient G is replaced by a smoothed value G in the descent process through a second-order implicit formulation [44]. This acts as a preconditioner which allows the use of much larger step lengths and ensures that each new shape in the optimization sequence remains smooth. 4. Results

The terms multiplied by the variation dwni;j;k of the discrete flow variables from Eq. (8) are collected and are equated to zero in order to eliminate them from the Lagrangian function to establish the discrete adjoint fluxes,

The first part of the results section contains a code validation of the NLFD approach for the solution of the Caradonna–Tung hover test cases. The study compares the pressure distributions at several

7

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15

−0.5 NLFD − 3 SPP Experimental Data [38]

Pressure Coefficient, Cp

span stations for three different test cases. In the first test case, the collective pitch was set at 0° and the tip Mach number at 0.520. For the second and third cases the tip Mach numbers are 0.439 and 0.877 at a collective pitch of 8°. Convergence history of the NLFD flow solver for all three test cases are illustrated as well. In the second part of this section, a redesign of the Caradonna–Tung blade at two different tip Mach numbers are demonstrated. The tip Mach numbers are set to 0.860 and 0.877 with a collective pitch of 8°. Initially the convergence of the NLFD adjoint solver is demonstrated and the gradient for various time steps are compared at three spanwise locations to establish the minimum number time steps required for design optimization. The initial and final pressure distributions, as well as the convergence of the objective function and blade cross-sections, are compared and discussed.

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x/c Fig. 6. Comparison of surface pressure distribution at z/R = 0.80 (collective pitch of 0° and tip Mach of 0.520).

Pressure Coefficient, C

p

−0.5

0

0.5

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x/c Fig. 7. Comparison of surface pressure distribution at z/R = 0.89 (collective pitch of 0° and tip Mach of 0.520).

0

10

−0.5

Zeroth Mode First Mode

−1

10

Pressure Coefficient, Cp

log (Normalized Max Residual)

0.5

1

4.1. Validation against the Caradonna–Tung experimental blade The first test case is of a non-lifting rotor. The tip Mach number is set to 0.520 at a collective pitch of 0°. Fig. 5 demonstrates the convergence of the NLFD flow solver for both the zeroth and first modes for a three time step per period solution on a 256  64  48 C–H grid, where the introduction of higher modes are due to the interpolation at the periodic boundary surfaces. The final solutions were acquired after a five-order magnitude drop in the residual requiring approximately 500 multigrid cycles. Both the zeroth and first modes have similar convergence rates. A higher convergence rate is observed within the first 100 multigrid cycles and the rate reduces thereafter. The Courant–Friedrichs–Levy (CFL) number had to be dropped to five for both the fine and coarse grids during the multigrid cycle, instead of the usual eight for inviscid three-dimensional steady and unsteady calculations. This is mainly due to the higher aspect ratio cells in the vicinity of the blade tip at the periodic boundaries. For convergence acceleration a 4 W multigrid cycle as well as implicit residual smoothing were utilized. At the far-field, the choice of boundary conditions either Riemann invariants or hover boundary conditions did not affect the solution. Figs. 6–8 compares the pressure distribution against the experimental data with excellent agreement at all three span locations. The simulation was not computed for a higher number of time steps per period since the three time step case produced sufficient results due to the steady nature of the solution. The computed thrust coefficient, Ct, is 0.00001, which certainly is within an acceptable range for a non-lifting rotor. Simulations were performed on a cluster based on single core 2.2 GHz AMD Opteron

0

−2

10

−3

10

−4

10

0

0.5

−5

10

−6

10

0

50

100

150

200

250

300

350

400

450

500

Multigrid Cycles Fig. 5. Convergence history of NLFD flow solver (collective pitch of 0° and tip Mach of 0.520).

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x/c Fig. 8. Comparison of surface pressure distribution at z/R = 0.96 (collective pitch of 0° and tip Mach of 0.520).

8

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15 0

10

−1.5

NLFD−3 0−Mode NLFD−5 0−Mode

−1

10

NLFD−5 1−Mode NLFD−7 0−Mode

−2

10

NLFD−7 1−Mode NLFD−7 2−Mode NLFD−7 3−Mode

−3

10

−4

10

−5

10

−1

p

NLFD−5 2−Mode

Pressure Coefficient, C

log (Normalized Max Residual)

NLFD−3 1−Mode

−0.5

0

0.5

1

−6

10

−7

10

1.5 0

100

200

300

400

500

600

700

800

900

0

1000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x/c

Multigrid Cycles

processor with 2 GB of memory per socket and communication between nodes are enabled by a switch based on the Infiniband Standard. Total simulation time was 30 min on 12 processors. The next two cases are of a lifting rotor at a collective pitch of 8 °. The first case is at a tip Mach number of 0.439. The simulation was acquired on a 256  64  48 C–H grid. Fig. 9 illustrates the convergence history of the NLFD flow solver for three different cases, each with an increasing number of modes. The flow solver converges six-orders of magnitude within 1000 multigrid cycles for the zeroth modes of the three, five, and seven time steps per period cases and the higher modes are excited for all three cases and drop four-orders of magnitude. The convergence rates are very similar when the number of modes are increased. The zeroth modes are all almost identical and the same can be said of the higher modes. Similar to the zero degree case, the convergence rate is higher for the first 100 cycles, which is a typical behavior for explicit codes. The CFL number was set at five and hover boundary conditions are employed at the far-field boundaries. In Figs. 10– 12, a validation of the surface pressure coefficient is presented for the three, five, seven, and nine time steps per period cases. The figure examines the pressure distributions at three separate span stations, 80%, 89%, and 96% and compares the results to the

Fig. 11. Comparison of surface pressure distribution at z/R = 0.89 (collective pitch of 8° and tip Mach of 0.439).

−1.5

−1

Pressure Coefficient, Cp

Fig. 9. Convergence history of NLFD flow solver (collective pitch of 8° and tip Mach of 0.439).

−0.5

0

0.5

1

1.5 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x/c Fig. 12. Comparison of surface pressure distribution at z/R = 0.96 (collective pitch of 8° and tip Mach of 0.439)

−3

x 10

−1.5 NLFD − 3 TSpP

14

NLFD − 5 TSpP

−1

NLFD − 9 TSpP Experimental Data [38]

−0.5

0 −1 −0.8

0.5

NLFD − 3 TSpP NLFD − 5 TSpP

−0.6

NLFD − 7 TSpP −0.4

1

0

0.1

0.2

12 NLFD − 3 TSpP (C = 0.00462) t

10

NLFD − 5 TSpP (C = 0.00462) t

NLFD − 7 TSpP (Ct = 0.00462)

8

NLFD − 9 TSpP (Ct = 0.00460)

Exp. Data [38] − Ct = 0.00459

6

NLFD − 9 TSpP 0.02

1.5

Thrust Coefficient, Ct

Pressure Coefficient, C

p

NLFD − 7 TSpP

0.3

0.04

0.4

0.06

0.08

0.5

4

0.1

0.6

0.7

0.8

0.9

x/c Fig. 10. Comparison of surface pressure distribution at z/R = 0.80 (collective pitch of 8° and tip Mach of 0.439).

0

50

100

150

200

250

300

350

400

450

500

Multigrid Cycles Fig. 13. Convergence history of the thrust coefficient (collective pitch of 8° and tip Mach of 0.439).

9

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15 0

log (Normalized Max Residual)

10

NLFD−3 0−Mode NLFD−3 1−Mode NLFD−5 0−Mode NLFD−5 1−Mode NLFD−5 2−Mode NLFD−7 0−Mode NLFD−7 1−Mode NLFD−7 2−Mode NLFD−7 3−Mode

−1

10

−2

10

−3

10

−4

10

−5

10

−6

10

−7

10

0

100

200

300

400

500

600

700

800

900

1000

Multigrid Cycles Fig. 15. Convergence history of NLFD flow solver (collective pitch of 8° and tip Mach of 0.877).

−1.5 NLFD − 3 TSpP NLFD − 5 TSpP

−1

NLFD − 7 TSpP

p

Pressure Coefficient, C

Caradonna–Tung data [40]. The three figures demonstrate that the three time step solution was sufficient to acquire a very good comparison against experimental data. In Fig. 10 an inset illustrates the comparison among the four cases for the upper surface pressure distribution in the vicinity of the leading edge, where an almost point-to-point match is observed. In Fig. 13, a convergence history of the thrust coefficient for the various cases is illustrated. The convergence history is almost identical among the four cases and global convergence is attained at approximately within 185 multigrid cycles, where the flow solver has achieved a three-order of magnitude reduction in the residual. The three, five, and seven time step cases converge to a thrust coefficient, Ct of 0.00462, whereas the nine time step solution converged to 0.00460. These values compare very well against the experimental value of 0.00459. The nine time step per period case certainly is closer to the experimental value and one could make the case that these many time steps are required; however, since the difference is small, we believe that as expected three time steps are sufficient for this steady test case. In Fig. 14 a plot of the spanwise sectional lift distribution for the various temporal cases demonstrates a very good comparison against experimental values except at the blade tip, where higher lift is predicted. The second of the two lifting rotor cases was computed at a tip Mach number of 0.877. Fig. 15 demonstrates the NLFD flow convergence history for three different cases. The flow solution was acquired after 1000 multigrid cycles with a six-order reduction in the residual of the zeroth and higher modes. All cases converged at approximately the same rate. The CFL as well as all other solver parameters are the same as the previous case. The convergence curve unlike the test case with a tip Mach number of 0.439 contains oscillations. Figs. 16–18 compare the surface pressure distribution for various temporal discretizations against the experimental data. In this test case, a region of supersonic flow is present from approximately the 80% blade span station onwards to the blade tip. A comparison against the experimental data shows differences near the vicinity of the upper surface shock wave, but the lower surface and pre- and post-shock regions are captured very well. The upper surface shock is slightly stronger than what is presented in the experimental data. This is mainly due to the fact that the simulations presented in this work are inviscid, and the lack of a boundary layer leads to an inaccurate prediction of the shock location. At the 80% spanwise location, the shock is weaker and the peak pressure is not captured as well when compared to the subsonic cases. We believe that this is primarily attributable

NLFD − 9 TSpP Experimental Data [38]

−0.5

0

0.5

1

1.5

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x/c Fig. 16. Comparison of surface pressure distribution at z/R = 0.80 (collective pitch of 8° and tip Mach of 0.877).

0.4

−1.5

0.35

Pressure Coefficient, C

p

−1

Lift Coefficient, Cl

0.3 0.25 0.2

NLFD − 3 TSpP NLFD − 5 TSpP NLFD − 7 TSpP NLFD − 9 TSpP Experimental Data [38]

0.15 0.1

−0.5

0

−1.2

NLFD − 3 TSpP NLFD − 5 TSpP NLFD − 7 TSpP NLFD − 9 TSpP

−1 −0.8

0.5

−0.6 −0.4

1 0.05

−0.2 0.35

0 0.4

1.5 0.5

0.6

0.7

0.8

0.9

1

r/R Fig. 14. Comparison of spanwise sectional lift distribution (collective pitch of 8° and tip Mach of 0.439).

0

0.1

0.2

0.4

0.3

0.4

0.45

0.5

0.6

0.7

0.8

0.9

x/c Fig. 17. Comparison of surface pressure distribution at z/R = 0.89 (collective pitch of 8° and tip Mach of 0.877).

10

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15 −4

−1.5

x 10

5

78

−0.5

0 NLFD − 3 TSpP NLFD − 5 TSpP NLFD − 7 TSpP NLFD − 9 TSpP

−0.8

0.5

−0.6 −0.4

1

1.5

79

80

81

82

83

0

Gradient

Pressure Coefficient, Cp

−1

−5

−0.2

0.38

0.4

0.2

0.3

0.42

0.44

0.46

−10 91

0

0.1

0.4

0.5

0.6

0.7

0.8

92

93

94

0.9 20

x/c

40

60

80

100

120

Design Variable

Fig. 18. Comparison of surface pressure distribution at z/R = 0.96 (collective pitch of 8° and tip Mach of 0.877).

Fig. 21. Comparison of gradients for various number of time steps at z/R = 0.89 (collective pitch of 8° and tip Mach of 0.877).

0

10

−4

−2

10

2 −4

10

1 76

78

80

82

−6

10

Gradient

log (Normalized Max Residual)

x 10

3

NLFD−3 0−Mode NLFD−3 1−Mode NLFD−5 0−Mode NLFD−5 1−Mode NLFD−5 2−Mode NLFD−7 0−Mode NLFD−7 1−Mode NLFD−7 2−Mode NLFD−7 3−Mode

−8

10

−10

10

−12

10

0

100

200

300

−1 −2 −3 −4

−14

10

0

400

500

600

700

800

900

1000

Multigrid Cycles

92

20

Fig. 19. Convergence history of NLFD adjoint solver (collective pitch of 8° and tip Mach of 0.877).

93

94

95

−5 40

60

80

100

120

Design Variable Fig. 22. Comparison of gradients for various number of time steps at z/R = 0.96 (collective pitch of 8° and tip Mach of 0.877).

−4

x 10

6 4

NLFD − 3 TSpP NLFD − 5 TSpP NLFD − 7 TSpP NLFD − 9 TSpP

Gradient

2 0 −4

x 10

−2

4 2

−4

0 −2 −4

−6

−6 80

85

90

95

−8 20

40

60

80

100

120

Design Variable Fig. 20. Comparison of gradients for various number of time steps at z/R = 0.80 (collective pitch of 8° and tip Mach of 0.877).

to the fact that the vortex wake is diffused in the far-field due to the large grid sizes and its induced effect on the 80% and inboard span stations are reduced. A grid independent study was not conducted in this work and is part of our future work. The solutions for the various temporal discretizations illustrate that the three time step per period is sufficient to acquire accurate pressure distributions. However, at both the 89% and 96% span stations, an inset shows a comparison of the four different cases at the location of the shock and a distinctive difference can be observed. Both the three and five time step cases show very similar results but the seven time step case establishes a slightly weaker shock. An additional simulation with nine time steps confirms the solution as it demonstrates an almost point-to-point match with the seven time step simulation at both spanwise locations. Therefore, for this particular transonic case, a minimum of seven time steps per period is required to obtain a temporal converged solution. The thrust coefficient converged to 0.00512, which is higher than the experimental value of 0.00473. This is primarily due to the higher sectional lift coefficients in the region of the blade tip.

11

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15 0

10 C /C T

To

C /C Q

−1

10

Qo

0.95

ΔI

Objective Function Ratio

1

−2

10

0.9 −3

10

0.85

−4

0

5

10

15

20

10

25

0

5

10

15

20

25

Design Cycle

Design Cycle

Fig. 23. Convergence of thrust and torque at a collective pitch of 8° and tip Mach of 0.877.

0

10

CT /CTo

0.98

−1

CQ/CQo

0.96

10

ΔI

Objective Function Ratio

1

0.94

−2

10

0.92 −3

0.9

10

0

2

4

6

8

10

12

14

16

0

2

4

6

8

10

12

14

16

Design Cycle

Design Cycle

Fig. 24. Convergence of thrust and torque at a collective pitch of 8° and tip Mach of 0.860.

4.2. Redesign of the Caradonna–Tung experimental blade The design process used in this work will change the shape of the rotor blade in order to minimize a chosen figure of merit, subject to the constrain that the thrust coefficient is maintained. The NLFD adjoint based design procedure, as illustrated in Fig. 4, requires the following steps: First, a set of multigrid cycles is used to drive the unsteady residual to a negligible value for all the modes used in the representation of the periodic solution. Second, the NLFD adjoint equations are driven down until four to five-orders of magnitude drop is observed. Third, an integral over the last period of the adjoint solution is used to form the gradient. This gradient is then smoothed using an implicit smoothing technique and the rotor shape is then modified in the direction of improvement as described in the previous section. Fourth, the internal grid is modified based on perturbations on the surface of the blade. The method modifies the grid points along each grid index line projecting from the surface. The arc length between the surface point and the far-field point along the grid line is first computed, then the grid point at each location along the grid line is attenuated proportional to the ratio of its arc length distance from the surface point and the total arc length between the surface and the far-field. The entire design process is repeated until the objective function converges and the gradient converges to acceptable levels The problems in this work typically required between 15 and 40 design cycles to reach the optimum.

The objective function is a weighted sum of the difference of the current and initial thrust coefficients and the torque coefficient as written below

I ¼ -1 C q þ

-2 2

ðC t  C tinitial Þ2 ;

where -1 and -2 are the weights for each objective function. The magnitude of the weights were set based on trial and error, and the final values will be presented and discussed in the following sections. The initial thrust coefficient, C tinitial , is set at the end of the initial flow solver cycle as the goal of the objective is to maintain the current level of thrust while reducing the torque coefficient. This section documents the approach for the redesign of the Caradonna–Tung experimental blade to effectively reduce the torque coefficient, while maintaining the thrust at every phase of the period. Here we establish the effect of the constraints on the final design by imposing weights via trial and error. These coefficients are maintained through the inclusion of these terms within the objective function as shown in the beginning of the section. During the initial computation to acquire the first solution, the thrust coefficients are saved in memory at each phase of the period and are subsequently used in the objective function as target values during the design cycles. An alternate procedure would be to prescribe desired target thrust values if they are known apriori. Its important to mention that a sensitivity analysis of the base code has been performed in the past [43,29].

12

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15

Fig. 25. Comparison of initial and final pressure distributions at various span stations at a collective pitch of 8° and tip Mach of 0.877.

All design runs employed surface mesh points as the design variables. There are 129 points on the surface of the blade at each of the 32 span stations, totaling to 4128 design variables. The gradients are implicitly smoothed in the chord-wise as well as the span-wise directions. For each design run, 500 multigrid cycles were employed for the flow and adjoint solvers to compute the initial solutions and only 35 were used for both solvers thereafter at each design cycle. This was sufficient to ensure global convergence of both the thrust and torque coefficients. Typical design runs required between 4 and 7 h of wall clock time on 12 processors. The convergence history of the adjoint solver is shown in Fig. 19 for three cases of increasing number of modes. All cases

experienced a reduction of at least eight-orders of magnitude within 1000 multigrid cycles. The convergence rate is independent of the number of modes for the first 500 cycles, where a five-order drop is seen, but the convergence rate reduces thereafter for the three and five time steps per period cases. This may be due to the nonlinear coupling effect between the modes. For design optimization, the accuracy of the gradient is of vital importance and the gradient of the base adjoint solver has been extensively compared to that obtained via the finite-difference and complex step approaches [43,29]. In this work, we investigate the effect of the number of modes on the accuracy of the gradient. In Figs. 20–22, this effect is investigated and compared at three different span

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15

13

Fig. 26. Comparison of initial and final pressure distributions at various span stations at a collective pitch of 8° and tip Mach of 0.860.

stations. At the 80%, the gradient compares well for all four cases. The inset shows a very small difference in the values of the gradient between the 70th and 100th design variable. In Figs. 21 and 22, the same is true, however a closer examination between the 77th and 84th as well as between the 90th and 95th points, shows that at least seven time steps is required to obtain converged gradient values. Since the difference is minor, we conclude that three is sufficient, and thus the subsequent design runs will only use three time steps per period. Figs. 23 and 24 illustrate the convergence of the thrust and torque ratios for two different Mach numbers, 0.860 and 0.877, respectively. The 0.877 tip Mach number case was the primary test case; however, it was observed and demonstrated later in this

paper that the shock was not completely eliminated. The 0.860 tip Mach number case was later added to demonstrate that almost complete elimination of the shock while maintaining the thrust coefficient is possible at a slightly lower tip Mach number. For both test cases, the collective pitch remains at 8°, the torque minimization contribution is present with a weight of unity, and the thrust coefficient is maintained with a weight of 200. In our initial runs, the weight for the thrust and torque coefficients were set to 5.0 and 1.0. This combination allowed the torque to be reduced within 15 design cycles, but maintaining the thrust coefficient proved to be a challenge. Upon closer examination, it was noticed that the sensitivities due to the torque coefficient were two-orders of magnitude higher, especially for the design variables close to the

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15

0.1 0.05 0 −0.05 −0.1 −0.15 −0.2

Initial Configuration Final Configuration

y/c

y/c

14

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.1 0.05 0 −0.05 −0.1 −0.15 −0.2

0.9

Initial Configuration Final Configuration

0.1

0.2

0.3

0.4

0.1 0.05 0 −0.05 −0.1 −0.15 −0.2

Initial Configuration Final Configuration

0.1

0.2

0.3

0.4

0.5

0.5

0.6

0.7

0.8

0.9

x/c

y/c

y/c

x/c

0.6

0.7

0.8

0.1 0.05 0 −0.05 −0.1 −0.15 −0.2

0.9

Initial Configuration Final Configuration

0.1

0.2

0.3

0.4

x/c

0.5

0.6

0.7

0.8

0.9

x/c

0.1 0.05 0 −0.05 −0.1 −0.15 −0.2

Initial Configuration Final Configuration

y/c

y/c

Fig. 27. Comparison of initial and final airfoil shapes at a collective pitch of 8° and tip Mach of 0.860.

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.1 0.05 0 −0.05 −0.1 −0.15 −0.2

0.9

Initial Configuration Final Configuration

0.1

0.2

0.3

0.4

0.1 0.05 0 −0.05 −0.1 −0.15 −0.2

Initial Configuration Final Configuration

0.1

0.2

0.3

0.4

0.5

0.5

0.6

0.7

0.8

0.9

x/c

y/c

y/c

x/c

0.6

0.7

0.8

0.9

x/c

0.1 0.05 0 −0.05 −0.1 −0.15 −0.2

Initial Configuration Final Configuration

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x/c

Fig. 28. Comparison of initial and final airfoil shapes at a collective pitch of 8° and tip Mach of 0.877.

location of the shock. Increasing the weight on the thrust coefficient to 200 allowed for a better maintenance of the thrust coefficient, but the final design required a larger number of design cycles. As shown from Fig. 23, the thrust is conserved to within 0.15% of the original value of 0.00512 and deviated not more than 0.2%, without significantly jeopardizing the reduction of the torque. As for the torque coefficient, it required approximately 20 design cycles and several additional cycles were performed to ensure that the minimum has been achieved. The design runs were typically terminated when the change in the torque coefficient dropped by three-orders of magnitude. At the 25th design cycle the torque coefficient reduced by 13% for the tip Mach number of 0.877. As for the case, where the tip Mach number is 0.860, the torque coefficient reduced by 8.6% while the final thrust coefficient was within 0.12% and deviating not more than 0.26% of the initial value of 0.00470. Figs. 25 and 26 demonstrate the initial and final pressure distributions at four span stations located at 80%, 89%, 92%, and 96% for both the 0.877 and 0.860 tip Mach numbers. The 92% span station was added to the previous three as it was observed that the shock is strongest at this location. Since the objective is to reduce the inviscid torque coefficient, it would be important to observe the change in the pressure distribution at this span station. In Fig. 25, at both the 92% and 96% span stations, an almost complete

elimination of the shock is observed, and the peak pressure increases by a small amount. At the 89% span station, the shock is severely weakened but it has been translated aft and the peak pressure increases slightly. This is primarily due to the need to maintain the thrust coefficient. In the case, where the tip Mach number is 0.860 as shown in Fig. 26, a slightly weaker shock is generated on the upper surface. Similar to the 0.877 tip Mach number case, at the 92% and 96% span stations the pressure discontinuity is completely eliminated. At the 89% span station, the shock is almost completely eliminated as opposed to the 0.877 case, where the remnants of a weak shock is still visible as it is also seen in the final pressure contour plots for both cases. These test cases conclude that the pressure discontinuity can be severely weakened but not eliminated and the thrust coefficient maintained for a single design point optimization of a hovering rotor blade. Figs. 27 and 28 examine the initial and final blade cross-sections at the 80%, 89%, 92%, and 96% span stations. There are three distinctive features of the new cross-section. First, the reduction of the upper surface curvature of the first 30% of the chord at all span stations. Second, an increase in the thickness as well as curvature between the 40% and 60% of the chord-wise locations. The combined effect of the two shape modifications is to introduce camber into a symmetrical airfoil to maintain the thrust coefficient while reducing the torque. Third, a small reduction of the leading edge radius at all four span stations produced higher suction peaks to assist in

C.A. Tatossian et al. / Computers & Fluids 51 (2011) 1–15

the maintenance of the initial thrust level. The effect of these three shape changes is visible in the final pressure distributions as shown in Figs. 25 and 26. 5. Conclusion An aerodynamic shape optimization approach has been developed to improve the performance of hovering rotor blades in transonic flow using a discrete adjoint method based on the Non-Linear Frequency Domain approach. The flow solution was validated against experimental data at both non-lifting at a tip Mach number of 0.520 and collective pitch of 0° and lifting cases at tip Mach numbers of 0.439 and 0.877 and collective pitch of 8°. As expected the NLFD technique produced sufficiently accurate results for both the adjoint and gradient vector with just three time steps per period. The discrete adjoint based NLFD framework proved to be beneficial in reducing the inviscid torque coefficient, while maintaining the initial thrust coefficient. At a tip Mach number of 0.860, the upper surface shock was almost completely eliminated, resulting in a 8.6% reduction in the torque coefficient while maintaining the thrust coefficient within 0.12% of the initial value. At a slightly higher tip Mach number of 0.877, the stronger upper surface shock was severely weakened, resulting in a reduction of 13% of the initial torque coefficient, while maintaining the thrust at 0.00512. Typical design run times were between 4 and 7 h on 12 processors. These results demonstrate a proof-of-concept of the potential of the approach to provide improvements to the performance of aerodynamic shapes that are subjected to unsteady loads. The extension of the approach for rotors in forward flight, viscous flows, and the possibility of including rotor planform parameters as design variables is part of our future work. Acknowledgment This research has benefited from the generous support of the Natural Sciences and Engineering Research Council of Canada from both an NSERC PGS scholarship as well as NSERC Discovery Grant # 206260. References [1] Jameson A. Computational aerodynamics for aircraft design. Science 1989;245:361–71. [2] Reuther J, Cliff S, Hicks R, van Dam CP. Practical design optimization of wing/ body configurations using the Euler equations. AIAA Paper 92-2633; 1992. [3] Gallman J, Reuther J, Pfeiffer N, Forrest W, Bernstorf D. Business jet wing design using aerodynamic shape optimization. In: 34th Aerospace sciences meeting and exhibit, Reno, Nevada; January, 1996. AIAA Paper 96-0554. [4] Reuther J, Alonso JJ, Martins JRRA, Smith SC. A coupled aero-structural optimization method for complete aircraft configurations. In: 37th Aerospace sciences meeting and exhibit, Reno, Nevada; January, 1999. AIAA Paper 990187. [5] Kasidit L, Jameson A. Case studies in aero-structural wing planform and section optimization. In: 22nd AIAA applied aerodynamics conference, providence, RI; August, 2004. AIAA Paper 2004-5372. [6] Nadarajah S, Jameson A, Alonso JJ. An adjoint method for the calculation of remote sensitivities in supersonic flow. In: 40th Aerospace sciences meething and exhibit, Reno (NV); January, 2002. AIAA Paper 2002-0261. [7] Mavriplis D. Multigrid solution of the discrete adjoint for optimization problems on unstructured meshes. AIAA J 2006;44:42–50. [8] Soto O, Lohner R, Yang C. An adjoint-based design methodology for CFD problems. Int J Numer Methods Heat Fluid Flow 2004;14(6):734–59. [9] Giles MB, Duta M, Muller J, Pierce N. Algorithm developments for discrete adjoint methods. AIAA J 2003;15(5):1131–45. [10] Nemec M, Zingg D, H Pulliam T. Multipoint and multi-objective aerodynamic shape optimization. AIAA J 2004;42(6):1057–65. [11] Mohammadi B, Pironneau O. Shape optimization in fluid mechanics. Ann Rev Fluid Mech 2005;36:255–79. [12] Srinivasan GR, Baeder JD, Obayashi S, McCroskey WJ. Flowfield of a lifting rotor in hover: a Navier-Stokes simulation. AIAA J 1992;30(10):2371–8. [13] Ahmad JU, Strawn RC. Hovering rotor and wake calculations with an oversetgrid Navier-Stokes solver. In: 55th Annual forum of the american helicopter society; 1999. [14] Datta Anubhav, Sitaraman J, Chopra I, Baeder JD. CFD/CSD prediction of rotor vibratory loads in high-speed flight. J Aircraft 2006;43(6):1698–709.

15

[15] Hariharan N, Sankar LN. High-order essentially non-oscillatory schemes for rotary-wing wake computations. J Aircraft 2004;41(2):258–67. [16] Yeshala N, Egolf TA, Vasilescu R, Sankar LN. Application of higher order spatially accurate schemes to rotors in hover. In: 24th Applied aerodynamics conference, AIAA; June, 2006. Number AIAA 2006-2818. [17] Pomin H, Wagner S. Navier-Stokes analysis of helicopter rotor aerodynamics in hover and forward flight. AIAA Paper 2001-0998, 2001. [18] Gordon Lieshman, Principles of helicopter aerodynamics, 2nd ed., Cambridge Aerospace Sciences; 2006. [19] Nadarajah S, Jameson A. Optimum shape design for unsteady flows with a time accurate continuous and discrete adjoint method. AIAA J 2007;45(7):1478–91. [20] Siva Nadarajah. The discrete adjoint approach to aerodynamic shape optimization. Ph.d. Dissertation, Stanford (CA): Department of Aeronautics and Astronautics, Stanford University; January, 2003. [21] Jameson A. Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings. In: AIAA 10th computational fluid dynamics conference, Honolulu (HI); June, 1991. AIAA paper 91-1596. [22] Lee SW, Kwon OJ. Aerodynamic shape optimization of hovering rotor blades in transonic flow using unstructured meshes. AIAA J 2006;44(8):1816–25. [23] Morris AM, Allen CB, Rendall TCS. Development of generic cfd-based aerodyanmic optimization tools for helicopter rotor blades. AIAA Paper 2007-3809; 2007. [24] Adamczyk JJ. Model equation for simulating flows in multistage turbomachinery. Technical report, NASA Technical Memorandum 86869; November, 1984. [25] Hall KC, Thomas JP, Clark WS. Computation of unsteady nonlinear flows in cascades using a harmonic balance technique. Technical report, In: 9th International symposium on unsteady aerodynamics, aeroacoustics and aeroelasticity of turbomachines, Lyon, France; September, 2000. [26] McMullen M, Jameson A, Alonso J. Acceleration of convergence to a periodic steady state in turbomachinery flows. In: 39th Aerospace sciences meething and exhibit, Reno (NV); January, 2001. AIAA Paper 01-0152. [27] McMullen M, Jameson A, Alonso J. Application of a non-linear frequency domain solver to the euler and navier-stokes equations. In: 40th Aerospace sciences meething and exhibit, Reno (NV); January, 2002. AIAA Paper 02-0120. [28] Matthew McMullen. The application of non-linear frequency domain methods to the euler and Navier-Stokes equations. Ph.d. Dissertation, Stanford (CA): Department of Aeronautics and Astronautics, Stanford University; March, 2003. [29] Nadarajah S, McMullen M, Jameson A. Optimum shape design for unsteady flows using time accurate and non-linear frequency domain methods. In: 16th Computational fluid dynamics conference, Orlando (FL); June, 2003. AIAA Paper 2003-3875. [30] Nadarajah S, McMullen M, Jameson A. Non-linear frequency domain method based optimum shape design for unsteady three-dimensional flows. In: AIAA 44th aerospace sciences meething and exhibit, Reno (NV); January 2006. AIAA Paper 2006-1052. [31] Nadarajah S, Jameson A. Optimum shape design for unsteady threedimensional viscous flows using a non-linear frequency domain method. AIAA J Aircraft 2007;44(5):1513–27. [32] Nadarajah S. Convergence studies of the time accurate and non-linear frequency domain methods for optimum shape design. Int J Comput Fluid Dynam 2007;21(5-6):189–207. [33] Choi S, Alonso J, Weide E, Sitaraman J. Validation study of aerodynamic analysis tools for design optimization of helicopter rotors. AIAA Paper 20073929; 2007. [34] Kumar M, Murthy VR. Efficient rotor flow calculations using multiblade coordinate transformed cfd in frequency domain. AIAA Paper 2008-987; 2008. [35] Butsuntorn N, Jameson A. Time spectral method for rotorcraft flow. AIAA Paper 2008-403; 2008. [36] Duta MC, Giles MB, Campobasso MS. The harmonic adjoint approach to unsteady turbomachinery design. Int J Numer Methods Fluids 2002;40:323–32. [37] Thomas JP, Hall KC, Dowell EH. A discrete adjoint approach for modelling unsteady aerodynamic design sensitivities. In: 41th Aerospace sciences meething and exhibit, Reno (NV); January, 2003. AIAA Paper 03-0041. [38] Jameson A. Analysis and design of numerical schemes for gas dynamics 1, artificial diffusion, upwind biasing, limiters and their effect on multigrid convergence. Int J Comput Fluid Dyn 1995;4:171–218. [39] Holmes DG, Tong SS. A three-dimensional euler solver for turbomachinery blade rows. J Eng Gas Turb Power 1985;107:258–64. [40] Caradonna FX, Tung C. Experimental and analytical studies of a model helicopter rotor in hover. In: Sixth European rotorcraft and powered lift aircraft forum, September 16–19, 1980. [41] Message Passing Interface Forum. MPI: a message-passing interface standard, May, 1994. . [42] Nadarajah S, Jameson A. A comparison of the continuous and discrete adjoint approach to automatic aerodynamic optimization. In: 38th Aerospace sciences meeting and exhibit, Reno (Nevada); January, 2000. AIAA paper 2000-0667. [43] Nadarajah S, Jameson A. Studies of the continuous and discrete adjoint approaches to viscous automatic aerodynamic shape optimization. In: 15th Computational fluid dynamics conference, Anaheim (CA); June, 2001. AIAA paper 2001-2530. [44] Jameson A, Vassberg JC. Studies of alternative numerical optimization methods applied to the brachistochrone problem. OptiCON 1999;99.