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

For both the 2D and 3D cases optimizations are carried out at Reynolds numbers of 2, 000. This Reynolds number was chosen to avoid the above-mentioned problems associated with higher Reynolds numbers while still being well within the range seen in natural ﬂiers. This Reynolds number physically corresponds to creatures on the scale of large insects and small birds such as hummingbirds. Modern unsteady, viscous ﬂow solvers are well-suited and validated for this Reynolds number regime.

Another physical parameter choice that is highly relevant to numerical simulation is the Mach number. Flapping ﬂight near ground level in Earth’s atmosphere for both natural ﬂiers and for any man-made craft inevitably be in the very low Mach number range on the order of M = 0.01 and can be considered e↵ectively incompressible.

The assumption of compressible vs. incompressible ﬂow is an important distinction for ﬂow solvers. Very low Mach numbers towards the limit of incompressibility leads to the decoupling of the pressure terms in the Navier-Stokes equations. There are a number of methods to address this decoupling, including low-Mach preconditioning and solving the decoupled incompressible Navier-Stokes equations. All of these methods come with a non-trivial increase in computational cost compared with solving the coupled equations in the compressible regime and can also be less robust when running many di↵erent cases without problem-speciﬁc tuning.

For these reasons, a compressible ﬂow solver is used for all cases considered herein, and the free-stream mach number is ﬁxed at M = 0.2. This Mach number was chosen because it is the smallest Mach number that can be reliably used with a compressible ﬂow solver without running into numerical issues related to the decoupling of the CHAPTER 2. METHODOLOGY 7 pressure terms. This decision was made largely for the beneﬁt of the speed and reliability of the compressible formulation and at the cost of some physical accuracy.

Flow at M = 0.2 is considered only mildly compressible, however there is still the potential for more pronounced compressible e↵ects when examining a ﬂapping wing.

Flow velocities can be accelerated signiﬁcantly around the leading edge of a ﬂapping wing if the wing itself is moving quickly relative to the free-stream ﬂow. With a freestream Mach number of M = 0.2 maximum ﬂow velocities could approach transonic speeds, whereas the maximum ﬂow velocity in a M = 0.01 free-stream would likely be only mildly compressible. Examination of the results presented later in the document indicates that in practice the maximum Mach number encounters does not exceed M = 0.5, as is illustrated in the visualization of a representative ﬂapping wing test case shown in ﬁgure 2.4.

2.1.1 Dimensional Analysis and Scaling Laws There are a number of useful non-dimensional parameters that apply to the freestream ﬂow, the kinematics of ﬂapping and the interaction between the ﬂuid and the airfoil/wing. These are the Reynolds number, the Strouhal number and the reduced frequency.

Reynolds Number

For the Reynolds number associated with the motion of the wing tip, the reference velocity is taken to be the mean tip speed Uref = Utip = !R where ! is the mean angular rate of the wing and R is the wing semi-span length. The angular rate is in turn deﬁned as 2⇥f where ⇥ is the angle swept from the top to bottom of the stroke and f is the oscillation frequency. The reference length for the wing tip Reynolds number is again taken to be cm. The tip Reynolds number is thus.

Reduced Frequency The reduced frequency k relates to the degree of wing-induced aerodynamic unsteadiness in a given ﬂow and is generically deﬁned as

2.2 Optimization Framework A block diagram of the software framework for solving ﬂapping-wing optimization problems is shown in ﬁgure 2.5. The framework is logically broken into the three parts shown, the optimization algorithm, the ﬂow solver and the wrapper routines that allow the ﬂow solver and the optimization algorithm to interface.

The optimization algorithm maintains a set of parameters that are varied throughout the optimization process to ﬁnd the maximum achievable value for the prescribed performance metric. These parameters encode the kinematics of the 2D oscillation or 3D ﬂapping motion. The wrapper takes these parameters and converts them to input compatible with the ﬂow solver. The ﬂow solver then computes the time-accurate ﬂow solution for the given motion, computing various performance metrics, such as lift, drag and power, at each time-step. The wrapper then takes these performance metric data and integrates them to arrive at the objective function value, which is passed back to the optimization algorithm. The details of these individual components are described in the remainder of this chapter.

CHAPTER 2. METHODOLOGY 13

** Figure 2.5: Software block diagramCHAPTER 2. METHODOLOGY 14**

2.3 2D and 3D Flow Solvers The 2D and 3D ﬂow solvers are based on the low-dissipation kinetic energy preserving (KEP) ﬁnite volume scheme developed by Jameson [13, 14] and extended by Allaneau and Jameson [2]. The kinetic energy preserving property of this scheme allows stability to be maintained with little or no artiﬁcial dissipation. This property is especially desirable for vortex dominated ﬂows such as ﬂapping ﬂight since artiﬁcial dissipation tends to quickly and unnaturally dampen complex ﬂow features. It is also straight-forward to implement and requires no additional computation overhead.

These codes have been speciﬁcally developed to simulate oscillating airfoils and ﬂapping wings and to be integrated into an optimization framework. The full details of these codes can be found in the work of Allaneau et. al.[1].

**2.3.1 Numerical FormulationKinetic Energy Preserving Finite Volume Scheme**

The ﬂow solvers used to compute the aerodynamic performance metrics that in turn are used to compute the objective functions used in the optimization process are based on the Navier-Stokes equations in 2D and 3D. The governing equations can be written in conservation form as

where the last summation term accounts for the boundary contributions.

Each interior face appears twice in the sums on the right-hand side of (2.19) with opposite signs for its discrete face area, and thus the contribution is

Summary of KEP Conditions To summarize, the discrete global variation law for the conservation of kinetic energy i

**can satisﬁed if the elements of fop satisfy the following conditions:**

Time Integration Method A TVD Runge-Kutta second-order multi-stage scheme is used for time-stepping in the 2D and 3D codes [30]. For a semi-discrete equation of the form

do not guarantee preservation of kinetic energy in time. The semi-implicit CrankNicholson scheme preserves kinetic energy in time, however the afore-mentioned computational costs would be too great to be practical in the context of optimization.

2.3.2 Mesh Motion Mesh motion for the 2D cases is applied in a straight-forward manner by rigid-body rotation and translation of the nodal coordinates such that the mesh is rotated to the desired pitch angle and then translated according to the desired degree of plunge.

Mesh motion for the 3D cases cannot be achieved in such a straight-forward manner because wing surface deformations are required. We have chosen to implement a computationally-inexpensive analytic mesh motion scheme. The base mesh in 3D is a structured H-C mesh with nodal coordinates x(i, j, k), y(i, j, k) and z(i, j, k) where i = 1..imax, j = 1..jmax, k = 1..kmax.

** Figure 2.6 shows the boundary cells of an example mesh.**

The H-C mesh is constructed by stacking C meshes along the Z-axis, resulting in XY-planes being parallel.

We allow a ne transformations that are functions of z only. This guarantees that XY-planes remain parallel, avoiding intersections that would result in collapsed and negative volume cells.

**The set of allow mesh transformations includes the following:**

2.4 2D and 3D Meshes The 2D ﬂow solver uses a structured C-mesh. It is possible for the 2D cases to resolve to the Kolmogorov scales at a Reynolds number of Re = 2, 000 using a mesh of 4096 ⇥ 512 cells. The ideal case from the standpoint of physical accuracy would be to use this mesh throughout the optimization process, but in practice this mesh requires on the order of 5 hours per ﬂapping cycle on 256 processors and thus would require a signiﬁcant amount of time within the optimization framework. We thus use a mesh with 1024 ⇥ 128 cells during the optimization process and then verify the resulting motion using the DNS mesh. Details of the meshes used are illustrated in ﬁgure 2.8.

The 3D ﬂow solver uses and H-C mesh that is essentially a series of 2D C-meshes stacked in parallel along the Z-axis. A lift and drag convergence study was undertaken to determine the dimensions of the 3D mesh to be used for the optimizations. Three simulations of a ﬂapping wing were run using meshes with 128 ⇥ 32 ⇥ 32, 256 ⇥ 64 ⇥ 64 and 384 ⇥ 96 ⇥ 96 cells. Plots of lift and drag versus time are shown in ﬁgure 2.9.

The results show that there is a relatively signiﬁcant change in moving from the 128 ⇥ 32 ⇥ 32 mesh to the 256 ⇥ 64 ⇥ 64, but relatively little change further increasing the mesh resolution to 384 ⇥ 64 ⇥ 64. The 256 ⇥ 64 ⇥ 64 mesh requires approximately one hour per ﬂapping cycle using 256 computing cores, which is adequate for the optimization process.

In the optimization process the concern is with time-averaged values of force coe cients and power integrated over a ﬂapping period. The table below shows average values of drag coe cient and the power in the x, y and z directions integrated over the second period.

** Figure 2.8: 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.

CHAPTER 2. METHODOLOGY 28

2.5 Convergence to Quasi-Periodicity Flow solutions are advanced in time starting from a steady-state solution. The motion is allowed to linearly accelerate into the periodic ﬂapping motion over the course of one-half period in order to avoid any discontinuous motion at t = 0 where t is the solver time. At the end of the half-period acceleration the solver reaches T = 0, where T is the period count. Force and power quantities are integrated over a single ﬂapping cycle, so it is important that these aerodynamic quantities have reached a quasi-steady state that consistently repeats from period to period. Several test cases were run to evaluate how many ﬂapping cycles are required before a quasi-steady state is reached. Figure 2.10 shows lift-versus-drag polars for one of the 2D cases. Figures

2.12 and 2.13 show the period-to-period convergence of the lift and drag coe cients in 2D and 3D, respectively. Using the lift coe cient as an example, the period-to-period CL (T ) CL (T 1) convergence is deﬁned as where T is the ﬂapping period number and CLmax CLmax is the highest measured CL value.

For the 2D case the plots indicate that the metrics become quasi-period to within 5% in drag and under 1% in lift by the start of the third period and the drag is further reduced to within 1% by the start of the ﬁfth period. For the 3D case the plots indicate that the drag is within 5% of quasi-period at the beginning of the second period, the lift is within 1% of quasi-periodic at the beginning of the second period, and the drag is further reduced to within 1% by the beginning of the ﬁfth period. These results thus indicate that function evaluations can be performed with as little as three ﬂapping cycles in 2D and two ﬂapping cycles in 3D.

CHAPTER 2. METHODOLOGY 31

2.6 Choice of Optimization Algorithm Each evaluation of the objective functions requires the computation of a ﬂow solution over a quasi-stable period of the ﬂapping cycle in order to integrate the appropriate force and power values. It is critical to use an optimization algorithm that minimizes the number of objective function evaluations in order to make any attempts at optimization possible in a reasonable amount of time due to the already signiﬁcant computational cost of a ﬂow solution.

The initial decision involved in selecting an optimization algorithm is whether to use a gradient-based or a gradient-free algorithm. Gradient-based algorithms inherently require fewer total function evaluations provided that the objective function is smooth in the vicinity of the optimal point. Gradient-based algorithms are especially e cient in terms of function calls if the initial point is close to the optimal point, so the choice of starting conditions is important. A drawback to gradient-based methods is that the algorithm only guarantees a local rather than a global minimum, which can cause problems for highly multi-modal objective functions. A work-around for this issue is to run multiple optimization from di↵erent initial conditions. In any case, the cost of gradient-free methods make them infeasible for the current application.

We have thus chosen to use a gradient-based optimization algorithm as described in the remainder of this section.

Gradient-based based algorithms, of course, require that a gradient be computable.