Abstract
Keywords
Introduction
MAV drones are revolutionizing the world with their versatility and controllability. They are employed in the civilian sector by film producers, photographers and, in the near future, by enterprises for delivering goods in urban areas. They share the airspace with helicopters and airplanes but the operability costs and the dimensions of the latter make MAVs a real asset in difficult environments and urban areas (Figure 1 shows an example of a professional drone).

The Parrot ANAFI Ai. Courtesy of Parrot Drones.
For this reason, MAV usage is steadily increasing every year and they need to comply with new and more restrictive international regulations that include noise emissions and safety. Since in the future MAVs will fly above people, their noise emission should respect different criteria. Because of the high rotation speed of MAV propellers, their noise is unpleasant for the human ear and this is mainly caused by the highly coupled aerodynamic and aeroacoustic interactions between the rotors1,2, and with the drone body 3 .
Another concern, this time for mission capabilities, comes from their low endurance. The viscous drag induced by the low Reynolds number at which MAV rotors operate reduces the aerodynamic efficiency of the rotors, hence reduces the endurance for a given energy storage.
The objective of this paper is to tackle both problems by providing an optimization framework that complies with industrial time and cost constraints. A Non-Linear Vortex Lattice Method, firstly introduced in sections Non-Linear Vortex Lattice Method and Tonal noise - Farassat Formulation-1A, then validated in section Validation, is used in an optimization loop that exploits the multi-objectives properties of the python framework Pymoo (see section Pymoo package) to design acoustically and aerodynamically optimized rotors. The optimization results are presented and analyzed in section Results.
Numerical methods
Non-Linear Vortex Lattice Method
The Non-Linear Vortex Lattice Method used in this work has been previously presented by Jo et al4,5 for rotor flow predictions. It is based on the incompressible
The rotor wake is also modeled using vortex rings following a prescribed wake geometry. The sectional linear lift is obtained by applying the Kutta-Joukovsky theorem (

On the left, an unsteady simulation snapshot showing the blade lattices, wake lattices and vortex particles. On the right, a steady simulation snapshot showing the blade lattices with the prescribed wake lattices.
Tonal noise - Farassat Formulation-1A
The Formulation-1A presented by Farassat
7
has been implemented to the aforementioned Non-Linear Vortex Lattice Method code to capture the tonal noise spectrum emitted by the rotor. The tonal noise, for low-Reynolds and low-Mach number rotors, is generated by two sources: The thickness noise generated by the displacement of the fluid due to the blade passage, which depends purely on the blade geometry (through the The loading noise that is dependent on the unsteady and steady pressure distributions on the blade (
The observer position plays an important role through the
Validation
To validate the code, unsteady simulations were compared with experiments from Gojon et al
8
obtained on a two-bladed, NACA0012-profiled, 10
In addition, the acoustic simulations were compared to the experimental results of Gojon et al
8
and Parisot-Dupuis et al
10
. As shown in Figure 3, differences lower than 3dB for microphones located below the rotor plane and between 3dB and 6dB for microphones located above the rotor plane were observed. Since MAVs are usually flown above people, acoustic simulation for a microphone at 30

On the left, the microphone locations are displayed; on the right, the comparison between the experimental data and the simulations is shown.
Optimization
The objective of this paper is to optimize the geometry of a MAV rotor in hovering conditions. The two optimization objectives are: The aerodynamic Figure-Of-Merit
11
Previous rotor aeroacoustic, multi-objectives optimizations12–17 have shown that, since these two objectives are in conflict, a single solution for this optimization that simultaneously improves both objectives does not exist. Instead, a series of optimal solutions that form a Pareto-front can be calculated. In this section, the optimization framework is presented (Pymoo package), followed by the optimization algorithm in Optimization algorithm, its parameters in Optimization setup and the design variables and constraints in Design variables, geometrical reconstruction and constraints.
Pymoo package
The Python-native optimization framework Pymoo 18 is used in this work. This framework allows to make both single- and multi-objectives optimizations, based on different algorithms, like GA, CMAES, NSGA-II, etc. The reader is referred to Blank and Deb 18 for further details.
Optimization algorithm
Genetic Algorithms (GA) are inspired by Charles Darwin’s theory of natural evolution and based on the survival of the fittest, but also on the appearance of crossover combinations and mutations that can lead to fitter successive generations. Since they work with a population of solutions, different set of solutions can be maintained throughout the optimization to reach global minima, while gradient-based optimizations can get stuck to local ones. They can also identify the set of Pareto optimal solutions. The main difference between GAs lies in the survival and selection methods.
The NSGA-II (Non-dominated Sorting Genetic Algorithm) 19 was chosen for this study. This algorithm is one of the first genetic algorithm and was opportunely modified from the first version of the code to improve the convergence speed and use the elitism as a way to increase the performance and prevent the loss of good solutions. The pseudo-algorithm is shown in Appendix NSGA-II pseudo-algorithm - Algorithm 1.
Optimization setup
As highlighted in the previous subsection, genetic algorithms like the NSGA-II used in this study require the number of candidates that will be part of the initial population. In this case, a population of 100 candidates has been chosen. The candidates are randomly generated in order to increase the diversity and the possibility of finding good candidates. The selection is made through the tournament selection method that selects a number of individuals, compares their fitness and selects the parents that will then be modified through the crossover operations and mutations to give birth to the new generation.
The crossover operation used is the ”Simulated Binary Crossover” (more details on Deb et al 20 ). This operation, also called recombination, combines the genetic data of different parents to create a newborn.
The mutation operation, instead, randomly modifies the parameters by taking into account a given probability. The mutation probability is here set to zero.
Design variables, geometrical reconstruction and constraints
The rotor geometry on which the optimizations are based is described in Table 1:
Non-optimized rotor parameters.
The design variables taken into account in the optimization are listed in Table 2 along with the minimum and maximum bounds.
Design variables used in this work and their minimum and maximum bounds.
A continuous function allows the optimization algorithm to select the values of each parameter within the minimum and maximum bounds (the number of possible combinations is set by machine precision). Once the ten variables are calculated by the optimization algorithm, the NVLM code builds the new geometry. The chord and twist distributions are defined as follow:
The root chord and pitch angle are respectively fixed at 0.025m and 10 The CP The chord length at the tip of the blade is defined by TIP
Once the interpolation function is defined by the previous points, the geometry is interpolated and 10 spanwise values of the twist angles and the chord lengths are calculated to generate the geometry.
A similar approach for the skew and rake (commonly known as winglet) distributions has been used:
Both skew and rake angles and their derivatives are equal to zero at the root of the blade; The TIP CP
Figure 4 shows, on top, examples of twist distributions interpolations. The upward pointing triangles represent the four control points with two different positions and two different values of the pitch angle, while the downward pointing triangles represent the tip value (for the sake of visibility it was kept constant in the plot). On the bottom examples of skew/rake distributions are presented: in this case both control points (represented by the upward pointing triangles) and angles at the tip (the downward pointing triangles) are changed and 4 different distributions are created.

Twist (on top) and skew (on the bottom) distribution definition from the control point (the upward pointing triangles) and the tip value (the downward pointing triangles).
The NACA0012 airfoil has been kept constant to limit the size of the parameter space.
To optimize both endurance and noise of a The algorithm selects the values of the design parameters and the An aerodynamic calculation at 4000RPM is run and the mean thrust is calculated; By making the assumption that the thrust follows an ideal quadratic function (Thrust A new aerodynamic calculation is run and, with it, an acoustic one in order to extract the two objective values FM and 1st BPF SPL. If both aerodynamic and aeroacoustic simulations are converged and the thrust is within acceptable bounds (
In the following section, the results of four different optimizations are presented: one optimization does not take into account the rake/skew angles (called ”Baseline” in Figure 5), one takes into account the rake distribution (but not the skew), while another one the skew (but not the rake), and one takes both into account. This helps in assessing the role of rake and skew on optimal solutions.

Four different Pareto-fronts showing the influence of the skew/rake design variables in the optimization loop.
Results
In this section the results from the aerodynamic and aeroacoustic optimizations of MAV rotors using the genetic algorithm NSGA-II and the non-linear vortex lattice method coupled with the Farassat Formulation-1A solution of the Ffowcs-Williams and Hawkings analogy are presented. As explained in the previous section, two calculations per iteration are needed, which means that for a total of 40 generations (with 100 candidates each), 8000 simulations need to be computed. For this reason, having a low computational cost per simulation is crucial and, therefore, the steady-simulation feature has been preferred over the unsteady counterpart. The NVLM main parameters are summed-up in Table 3.
Simulation parameters used for the validation of the NVLM code.
Figure 5 shows the Pareto-fronts obtained with four different optimization setups, including or not rake and skew distributions in the design variables.
The “Baseline” Pareto-front represents an optimization with only 6 design variables, since it does not include skew and rake variables into the optimization loop. “Skew” and “Rake” Pareto-fronts are 8-variables optimizations: they include the skew and rake distributions in the rotor design respectively. This allows to slightly improve the performance of the Pareto set of optimal rotors. It is shown that including the skew variables helps to move the Pareto-front towards better aerodynamic efficiencies but does not yield significant changes in acoustic performance. When the rake variables are introduced in the loop, instead, the optimization seems more appropriate to find aerodynamically ideal rotors and acoustically stealthier ones with 5dB difference with respect to the Baseline curve. The last optimization, depicted with the solid line, includes both rake and skew distributions into the design parameters and has the widest Pareto-front, with both aerodynamic and acoustic performance overtaking the ones given by the dashed line rake-included Pareto-front. For this reason, only this last optimization run will be presented and analyzed in more details.
The dots of Figure 6 scatter plot represent all the candidates of the 10-variables optimization that include both rake and skew distributions. The dark-grey ones, that represent the last generation candidates, are all very close to the Pareto-front and show that the optimization algorithm has converged to the Pareto set of optimal solutions. The figure also shows two linear trends: the dashed line (referred to as “1st Pareto line”) and the dash-dotted one (“2nd Pareto line”).

Scatter plot showing all the simulated candidates. The last generation candidates, the best aerodynamic efficient and the acoustically stealthier candidates (whose geometries are shown in Figure 12) are highlighted with different colours.
As described in subsection Tonal noise - Farassat Formulation-1A, the total tonal noise depends on two sources (the loading noise and the thickness noise). For this reason it is important to analyze the interaction between them to understand the previous graph. Figure 7 shows four differently coloured Pareto-fronts.

Four scatter plots showing the candidates coloured by different variables.
The main results are:
The rotation speed (Figure 7 - top-right scatter plot) changes along the Pareto-front. By following the 1st Pareto line from right to left, the rotation speed increases until the junction between the two Pareto lines. Here, the rotation speed drops suddenly, to increase again along the 2nd Pareto line. While the influence of a reduced rotation speed on the reduction of tonal noise is explicit from equations 3 and 2, these results indicate that it may be counterbalanced by other terms such as rotor solidity (shown in Figure 10 and defined as follows: The thickness noise (Figure 7 - bottom-right scatter plot) follows the same trend of the rotation speed. A steady increase from right to left of the first Pareto line is followed by a sudden decrease in the thickness noise level and, again, an increase along the second Pareto line; At a far-field microphone location, the loading noise (Figure 7 - bottom-left scatter plot) depends mainly on the magnitude of thrust force. However, even for geometries having the same thrust, differences in the order of 3dB are depicted in the graph and depend on the loading distribution, but also on the rotation speed of the rotor (see equation 3); The geometries colored by the 2nd BPF SPL (Figure 7 - top-left scatter plot) show that optimizing the 1st BPF does not affect, in the same way, the 2nd BPF.

Scatter plot showing the candidates colored by rotor solidity.
Overall, it can be understood from these results that while loading noise dominates the acoustic footprint on the first Pareto line, it is competed by thickness noise on the second Pareto line. An increase in the rotation speed can be compensated for by an increase in rotor solidity to achieve a given target thrust. Furthermore, an increase in rotor solidity increases thickness noise through its explicit contribution in equation 2. Because the thickness noise is weak compared to the loading noise on the first Pareto line, an increase in rotor solidity that significantly impacts the SPL and a reduction in SPL can be directly correlated with a decrease in RPM. Conversely, because the thickness noise is of the same order of magnitude than loading noise on the second Pareto line, an increase in rotor solidity contributes to an increase in SPL (see Figure 10). Hence SPL is not solely correlated with RPM, which explains why rotor geometries with lower acoustic footprint may not be obtained at minimum rotation speeds. However, it is stressed that reducing the 1st BPF peak does not necessarily comes with a reduction in the 2nd BPF peak.
Figures 8 and 9 show the chord lengths, as well as the pitch, rake and skew angles of the different Control Points (CP) and Tip (TIP). Here are additional insights on the Pareto optimal solutions:
The chord at the control point reaches the maximum value allowed to the algorithm. This maximum value is located near the root (see top-left scatter plot of Figure 11); The chord at the tip (Figure 8 - top-right scatter plot) is the key parameter to understand the relative trends of the two Pareto lines: the geometries with higher Figure-Of-Merit have a smaller chord at the tip, which contributes to reduce tip vortex strength. The right Pareto line, instead, presents larger chord lengths at the tip to increase rotor solidity, which may affect both rotation speed and enhance thickness/loading noise cancellation mechanisms; Rotors with larger aerodynamic efficiency have a higher pitch angle at the control point (see Figure 8 - bottom-left scatter plot) because they need to recover the thrust force that is lost by the smaller chord at the tip. The control points are located near the root (see Figure 11); The pitch at the tip (Figure 8 - bottom-right scatter plot) is almost constant everywhere, except on the left side of the first Pareto line. Here, in fact, a higher rotation speed compensates for the lower twist at the tip; The rake value is almost constant for all the geometries of the Pareto-front. It points downward, reaches the maximum value allowed in the optimization and is located very close to the tip. This probably has a double effect: it helps to reduce the torque induced by the tip vortex, and add an outward component to the loading vectors near the tip. The changed loading vectors modify the phase of the loading noise emitted by the rotor compared to a rotor without rake. Although it is not shown here for the sake of conciseness, it is seen that the phase change induced by the outward force component, modified the interaction with the thickness noise, hence there is a destructive interference. As a consequence, the total noise is reduced; The Pareto-front shows two different and opposite skew angles, but the sudden change does not happen on the junction between the two Pareto lines, it happens on the first Pareto line.

Four scatter plots showing the candidates coloured by design variable values.

Two scatter plots showing the candidates coloured by rake angle (on the left) and skew angle (on the right) at the tip.

Scatter plot showing the candidates coloured by the control point variables.
This optimization run shows how three design parameters are important to the improvement of the performance of MAV rotors. The chord length at the tip is the key parameter to understand the relative trends of the two Pareto lines. A smaller chord at the tip can lead to aerodynamically more efficient rotors. Larger ones to acoustically stealthier rotors. The rake and skew distributions play important roles in both aerodynamic and acoustic performance. The former, inducing a downward-pointing blade tip, reaches the maximum values allowed to the algorithm. The skew, instead, deforming the blade in the forward direction for acoustically stealthier rotors and backward for aerodynamically more efficient rotors. The two best geometries (whose performance are highlighted in Figure 6) are shown in Figure 12.

On the left: the rotor with highest Figure-Of-Merit, on the right: with the lowest acoustic footprint.
Conclusion
The approach presented in this paper uses a non-linear vortex lattice method, coupled with Farassat’s aeroacoustic tonal noise model and the genetic algorithm NSGA-II of Pymoo package to optimize the rotor geometries and find aerodynamically efficient and aeroacoustically stealth MAV rotors. The two optimization objectives chosen for this study are the Figure-Of-Merit describing the rotor aerodynamic efficiency and the BPF SPL peak for a microphone located at a far-field distance of 1.62m, at an angle below the rotor plane of 30
All the optimization runs take into account the twist and chord distributions in the generation of new geometries (the airfoil chosen is the NACA0012), for a total of six design variables. Two optimizations with two additional variables, namely the rake and skew, are independently added. A last optimization run, including all ten design variables, has been performed and has further improved the set of optimal solutions.
The combined effects of both rake and skew on aerodynamic and aeroacoustic performance allowed to reach efficient and stealth rotors. The Pareto-front presents two linear trends with different slopes, referred to as Pareto lines. The parameter that splits the two lines is the chord length at the tip. A smaller chord, in fact, is preferable for aerodynamically more efficient rotors while a larger chord yields stealthier rotors.
The negative rake concentrated at the tip pushes the whole Pareto-front to aerodynamically more efficient rotors. Coupling the skew modifications to the rake ones, creates rotors with even higher Figure-Of-Merit and allows to go further on the acoustic counterpart too.
The BPF SPL is only a part of the acoustic spectrum. The acoustic improvements here obtained do not imply improvements over the entire frequency spectrum. In fact, the 2nd BPF peak of the best acoustic configuration is not the lowest value obtained. In addition, if the sensitivity of the human ear is to be considered, the A-ponderation should be applied.
These results are obtained by using the steady simulation capabilities of the NVLM code. However, to validate all the presented results, experimental tests are being prepared in the ISAE-SUPAERO anechoic chamber. It is also worth noting that this optimization does not take into account the inertia of the blades, the rotation speed at which the motor is more efficient, and the rotor mass. These parameters must be taken into account into future optimizations of commercial MAV rotors.
Copyright
Copyright © Parrot Drones SAS 2021
