# «HIGH FIDELITY OPTIMIZATION OF FLAPPING AIRFOILS AND WINGS A DISSERTATION SUBMITTED TO THE DEPARTMENT OF AERONAUTICS & ASTRONAUTICS AND THE COMMITTEE ...»

With CFD it is possible to numerically compute gradients using an adjoint-based method [12, 6]. Adjoint-based methods have the advantage that the cost of a gradient computation is independent of the number of optimization parameters, leading to very e cient solutions for high-dimensional problems. Adjoint-based optimization methods are widely used for steady-state problems such as aerodynamic shape optimization to minimize transonic drag on airfoils and wings [12]. Adjoint-methods have met with more limited success for unsteady problems, however [26]. The computation of an unsteady adjoint requires that the unsteady ﬂow equations ﬁrst be integrated fully forward in time. The adjoint equations are then solved by integrating the adjoint equations backward in time using the previously-solved-for ﬂuid state values.

** CHAPTER 2. METHODOLOGY 36 The issue is that the adjoint equations require the state values at each time-step from the forward integration of the governing equations.**

For a 3D simulation on a mesh of size M ⇥ N ⇥ P and with T time-steps this amounts to 5 ⇥ M ⇥ N ⇥ P ⇥ T doubles-precision ﬂoating point values. Typical 3D solutions using an explicit solver and modestly-sized 3D mesh involve ⇡ 105 106 cells and 105 107 time-steps, resulting in a storage requirement on the order of 101 1 to 101 4 bytes. This is far more data than can be stored in RAM, so the solution must be stored on disk, resulting in a signiﬁcant increase of computational cost related to storing and reading this data.

The other option is to compute the gradients numerically via ﬁnite-di↵erences.

Finite-di↵erences require on the order of n additional function evaluations for an optimization problem with n parameters. The number of parameters in the optimization problem is thus a signiﬁcant factor in the total computational cost. In this work it has been found that a su ciently rich range of ﬂapping motions can be deﬁned with a relatively small set of parameters as described in sections 2.8 and 2.9, with the result that the ﬁnite di↵erence approach is computationally feasible. The objective function evaluations required by ﬁnite-di↵erence calculations are also independent of each other and so can be parallelized across the number of parameters, allowing further speedup given su cient computational resources.

CHAPTER 2. METHODOLOGY 37

2.7 SNOPT Optimization Algorithm We have chosen to use the SNOPT package [8] for all of the optimizations presented here. SNOPT is a widely used software library that uses a gradient-based sequential quadratic programming (SQP) algorithm and is designed for use on constrained, non-linear optimization problems. The SQP method involves sequences of major and minor iterations. The major iterations sequentially drive a set of parameters towards a point that satisﬁes the conditions for optimality and the constraints on the problem.

The update to the parameters in each major iterations comes from the solution of a quadratic programming (QP) sub-problem that is itself solved iteratively via a sequence of minor iterations.

Consider a constrained, non-linear optimization problem of the form

2.8 2D Parameterizations In the 2D case we consider single-mode sinusoidal pitching and plunging motions. In the general case these motions take the form

These motions are achieved by time-dependent rigid body transformations of the computational mesh.

The parameterization used in all cases considered herein consists of four optimization variables: the frequency f, the pitch and plunge amplitudes ↵ and h and the phase di↵erence between pitch and plunge. The base angle of attack ↵0 is set to zero in all cases, yielding motions that are symmetric about the axis of free-stream ﬂow.

2.9 3D Parameterizations The 3D motion parameterization is intended to mimic typical ﬂapping motions with combinations of varying dihedral angle at the wing root and twisting along the span.

The 3D wing is allowed three principle types of motion, described using standard aviation terminology: twisting, dihedral and sweep. Twisting corresponds to spanwise rotation about the Z-axis, dihedral corresponds to rotation about the X-axis and sweep corresponds to rotation about the Y -axis. All motions are symmetric about the XY -axis.

The motion of the wing is parameterized by using a wing skeleton with one or more control points. The skeleton forms a kinetic chain with each control point holding a set of angles for twist, dihedral and sweep. Denote this by Ci = {↵i, ✓i, ⇤i }, where Ci denotes the ith control point, ↵ is the twist angle, ✓ is the dihedral angle and ⇤ is the dihedral angle. The angles in Ci are relative to those in Ci 1, so that a value of zero maintains the angle from the prior control point. The wing skeleton is illustrated in CHAPTER 2. METHODOLOGY 40 ﬁgure 2.14 Smooth variation of the twist angle ↵(s), dihedral angle ✓(s) and the sweep angle ⇤(s) (where s is the span-wise arc-length parameter) are generated by ﬁtting piecewise cubic splines through the control points. This results in a set of shears and rotation as a function of Z that are used to deform the mesh in the method described in section 2.3.2. Examples of the deformed wing surface are shown in ﬁgures 2.15,

2.16 and 2.17.

The time-dependent ﬂapping motion is parameterized by prescribing sinusoidal variation to a set of N control points Ci. The equations of motion are of the form

where, and ' are phase parameters and f is the frequency of oscillation.

The total number of parameters for a wing with N control points is thus 6N + 1.

Sometimes it is desirable to use fewer control points for certain degrees of freedom to reduce the number of total optimization parameters. An example would be a case with 2 twisting control points, 1 dihedral and 1 sweep control point. This can be achieved by hard-coding some of the control point parameters to a ﬁxed value. An example hard-coding for the 2-1-1 twist-dihedral-sweep example would thus map to CHAPTER 2. METHODOLOGY 41 ✓3 ✓2 ✓1 (a) Span-wise dihedral and sweep kinematic skeleton joint type

2.10 Components of the Objective Functions and Constraints Objective functions are based on time-averages of various integrated force and power coe cients calculated by the ﬂow solver. These include the average lift CL, average thrust CT = CD and the aerodynamic power averages Px, Py and Pz. We also di↵erentiate between the thrust-producing power Px and the non-thrust-producing powers Py and Pz. In all cases these averages are computed by integrating over a single ﬂapping cycle after a suitable number of periods have elapsed to allow the ﬂow to reach a quasi- periodic state.

Propulsive e ciency is a commonly used performance metric in propellers and ﬂapping wing systems and the maximization of propulsive e ciency is an attractive objective function for the optimization of these systems. The propulsive e ciency is the ratio of thrust-producing power to the mechanical power required to ﬂap the wing or to pitch and plunge the airfoil. However, as noted by Jones and Platzer[17], the propulsive e ciency is a discontinuous function, as can be seen in ﬁgure 2.18.

Furthermore, the discontinuity occurs in the vicinity of the maximum propulsive efciency, making it a very di cult problem for many optimization algorithms. This discontinuity occurs as the wing transitions from energy production to energy extraction, i.e. ﬂutter. In the exact balance of zero power production no power is required to ﬂap the wing and the denominator of the standard propulsive e ciency formulation becomes zero. For this reason we use a “modiﬁed propulsive e ciency” ⌘m given by Px ⌘m = q (2.56) P x + Py + P z This formulation removes the problem of a discontinuous objective function since, for q all practical cases, Px + Py + Pz 0.

Another interesting optimization objective is the problem of constrained power minimization, which is relevant to steady, level ﬂight of a ﬂapping vehicle. In steady, level ﬂight, the lift is equal to vehicle weight and the thrust is equal to the vehicle drag. The minimization of power consumed in producing thrust in an important CHAPTER 2. METHODOLOGY 45 objective for small MAV-scale aircraft due to their limited energy storage capacity.

This can be cast as a constrained optimization problem of the form

** Figure 2.18: Discontinuous behavior of propulsive e ciency Chapter 3 Pitching & Plunging Airfoil Optimizations The ﬁrst cases considered are the unconstrained maximization of propulsive e ciency and the thrust-constrained minimization of input power, all for the case of a pitching and plunging airfoil.**

The pitching and plunging airfoil case is parameterized by frequency f, pitching amplitude ↵, plunging amplitude h and the phase di↵erence between pitching and plunging ✓.

Optimization of both pure plunging and combined pitching and plunging are considered for the maximization of propulsive e ciency. This allows for analysis of the beneﬁts of allowing the additional degrees of freedom of pitch and phase di↵erence.

**The parameter mappings for the maximization of propulsive e ciency are as follows:**

The pitching and plunging case is the only type of motion considered for the thrust-constrained power minimization optimizations. Power is minimized for a range of prescribed target thrust coe cient values and at a range of angles of attack. This allows for the identiﬁcation of trends in the motion parameters and the minimum

## CHAPTER 3. PITCHING & PLUNGING AIRFOIL OPTIMIZATIONS 48

**power required to achieve a given value of thrust.**

3.1 Maximization of Propulsive E ciency For the 2D pitching and plunging airfoil case we consider the maximization of the modiﬁed propulsive e ciency for a NACA0012 airfoil at a Reynolds number of 1850 based on the chord and the free stream velocity, and at a mach number of 0.2. The numerical solutions for all 2D cases use a 1024⇥128 C-mesh and are advanced through ﬁve periods to ensure that the force and power reach a periodic state. The force and power coe cients are integrated over the ﬁnal oscillation cycle to obtain the averaged quantities for computation of the objective function. The solver is run in parallel on a large cluster, typically using 64 compute cores per ﬂow solution. Compute times are on the order of 2 hours per ﬂow solution. The optimal case is re-run using a 4096⇥512 mesh for validation and further analysis of the optimal case. This mesh is e↵ectively of DNS resolution for the near-ﬁeld around the airfoil for the given Reynolds number.

Aspects of the 2D mesh are shown in ﬁgure 3.1.

The parameter values obtained for the maximization of propulsive e ciency are

**as follows:**

Case Frequency ↵ h ⌘m Plunging 2.63Hz 0.252 12.1% Pitching/Plunging 4.61Hz 19.78 0.212 67.03 31.4% The results from the optimizations show that an increase in modiﬁed propulsive e ciency from 12.1% to a maximum of ⌘m = 31.4% is achieved by adding the pitching and phase degrees of freedom. It is also interesting to note that the frequency nearly doubles for the pitching and plunging case, while the total plunge amplitude has decreased by roughly 15%. Figure 3.2 shows snapshots from the trajectories of the plunging and the pitching/plunging cases. Note that these plots have been normalized to one period, so the the di↵erence in frequency is not illustrated. The trajectories clearly show that the pitching motion serves to angle the airfoil towards the tangent of the trajectory. The e↵ect of this is to reduce the e↵ective angle of attack during

## CHAPTER 3. PITCHING & PLUNGING AIRFOIL OPTIMIZATIONS 49

the middle portion of the up- and down-stroke when compared to the pure plunging case.The lift versus drag polars in ﬁgure 3.3 give insight into the di↵erences in force production between these two cases. Non-lifting cases with ⌘m = 100% must all lie along the line CL = 0 since Py and Pz must be everywhere zero. The e ciency of a given ﬂapping cycle can then generally be inferred by observing the aspect ratio of the polar. The wider and the polar along CD and the narrower along CL, the greater the magnitude of ⌘m, and additionally, for ⌘m to be positive the polar must show a net thrust. This is seen in ﬁgure 3.3 where the pitching and plunging case shows a signiﬁcant increase in thrust production compared with the plunging case with only a relatively small increase in lift through the cycle.

** Figures 3.4 and 3.**

5 show the sequence of objective values computed by SNOPT during each of the optimizations. Note that each optimization begins with a period of pseudo-random function calls to establish an initial basis for the optimization. The plunging case reaches the neighborhood of the optimum within around 30 function evaluations and proceeds to compute ﬁnite di↵erence gradients to reﬁne the parameters. The pitching and plunging case, with four parameters, requires a greater number of function calls, and shows a more obvious trajectory towards the optimum. Here SNOPT makes several initial large improvements in the 20 to 40 function evaluation range, followed by steady reﬁnement towards the optimum. Note that the large jumps and apparent discontinuities occur during the line search phase as the step length is increased to check for the possibility of further improvement.

## CHAPTER 3. PITCHING & PLUNGING AIRFOIL OPTIMIZATIONS 50

Figure 3.1: 2D Meshes used for optimization and analysis. The left image shows the boundaries and extent of the C-Mesh. The center ﬁgure shows a portion of the 1024⇥128 mesh used during function evaluations by the optimization algorithm. The left ﬁgure shows a portion of the 4096 ⇥ 512 DNS mesh used for validation, analysis and ﬂow visualization.** Figure 3.2: Comparison of the relative motion of plunging (top) versus pitching and plunging (bottom).**

The individual frames of the motion are taken at ten equally spaced intervals in a single period of ﬂapping. Note then that the horizontal axes does not represent either the explicit time or space dimensions.