2.1. Aerodynamic Solver: Non-Linear Vortex Lattice Method
The NVLM code developed by Jo et al. [
19,
73] is used here to calculate the performance and the airflow properties induced by MAV rotors in hovering and isolated conditions (
Figure 1). Two different solvers of the same NVLM code are used: a steady one, which models the wake as an uncontracted helicoid, and an unsteady one, which models the wake by means of free vortex particles (
Figure 1). The former is used during the optimization phase because of its computational speed, while the latter is used for ad hoc verification of the performance of the optimal geometries.
Both solvers are based on the potential flow theory which assumes that the flow is incompressible
, inviscid
, and irrotational
. The last hypothesis leads to a conservative vector field, that is, the velocity vector
is equal to the gradient of a velocity potential, noted as
:
Equation (
1) satisfies Laplace’s equation
. Since Laplace’s operator
is linear, multiple velocity potentials can be linearly added to create more complex flows.
In the case of the VLM, the final flow-field is a combination of elementary vortex rings (or lattices) with constant strength
. Since the rotor blade is considered to be thin (with respect to the chord), the 3D blade camber line is discretized into
lattices (where
and
are the number of lattices along the chord and span, respectively) and the circulations of the vortex rings are noted as
[
74] (see
Figure 2). These vortices generate an induced velocity at all points in space and must satisfy the non-penetrating condition at the blade surface. As shown in
Figure 2, each vortex lattice follows the
–
rule [
74]: the leading vortex ring segment is placed at
of the geometric cell chordwise length, and a collocation point is placed at
of the chordwise cell length. In this point, noted
, the non-penetrating condition is applied. The velocity induced by a vortex ring at the collocation point
is computed by considering the contribution of all the four vortex edges thanks to the Biot–Savart law [
74].
The total number of lattices
N is equal to the number of rotors multiplied by the number of rotor blades and the number of lattices per blade
. As the circulations
are the unknowns of the problem, the influence coefficients
are initially computed by considering the velocity induced by the four vortex segments of lattice
q onto the collocation point of lattice
p, with
. For each lattice
p, the non-penetrating flow condition can be written as follows:
By moving the terms
to the right-hand side, the following system is obtained:
is the velocity induced by all wake vortices at the collocation point of lattice of index p, the velocity induced at lattice p by the rotation of the blade, is the freestream velocity vector, and the vector normal to the lattice.
A sectional non linear correction on
is then performed using the XFOIL [
75] lookup tables to account for the low Reynolds number induced non-linearities (for more information, please refer to Li Volsi [
76]).
As the wake affects the velocity induced on the blade, and therefore the rotor performance through the term , several methods have been developed to propagate the wake behind the trailing edge of the blades. As previously mentioned, two different wake models are used in this work: a steady one, which models the wake as a prescribed and uncontracted helicoidal wake, and an unsteady one, more accurate and computationally more expensive, that uses a vortex particle method (VPM).
On the latter, Lagrangian-described vortex particles (or
blobs) are shed at each timestep from the trailing edge of the blades. Their total number
depends on the number of blade spanwise lattices, the angle-defined timestep (referred to as
step angle afterwards), and the number of total revolutions. They are characterized by a constant circulation
and allow us to reconstruct the vorticity field at a given location
x anywhere near the rotor through the following formula:
with:
Since the
vortex particles constantly influence each others and their number increases throughout the simulation, the computational cost is proportional to
, but the implementation of a Fast Multipole Method has reduced it to
[
19].
The steady solver, on the other hand, is based on a prescribed wake model [
77]. The discretization of the helical wake is fixed from the beginning of the simulation and depends on the blade spanwise lattices, the
step angle (that, in this case, defines the angular spacing between two consecutive rows of wake lattices), and the number of revolutions to simulate. The induced velocity is calculated at each spanwise section
j by using Biot–Savart’s law and the blades’ lattices circulations. In the present prescribed wake model, the helix pitch is supposed to be constant from the inboard to the outboard part of the helix, the rotor induced velocity is calculated by averaging all the spanwise induced velocities and used to set the helix pitch. To satisfy the 3D Kutta condition at the trailing edge of the blades, the circulation of the first row (noted
k) of wake lattices must be equal to the circulations of the trailing edge lattices:
The other wake lattice rows take the same circulation strength as the first row:
is the total number of wake lattice rows.
In the next iteration, the new averaged spanwise-induced velocity projected on the rotational axis is calculated by taking into account the new blade and wake lattice circulations, so that the new helix pitch is calculated. As this is a steady simulation and the helix discretization is fixed at the beginning of the computation, the calculation stops when the thrust coefficient is converged.
The main outputs of both solvers are the rotor thrust and torque and the blade loading distribution (this is one of the inputs of the acoustic propagator presented in
Section 2.2). The aerodynamic efficiency is calculated through the figure of merit [
78] (FM, Equation (
7)), which depends on the thrust
T, torque
Q, rotational speed
, air density
, and rotor radius
R:
The simulation parameters, shown in
Table 1 and chosen for the optimizations, are sufficient to capture with a good accuracy both aerodynamic and aeroacoustic performance of the rotors, while limiting the computational time (please refer to Jo et al. [
19] and Li Volsi [
76] for the validation of the code). As previously explained, the
step angle defines the geometric discretization of the helical wake in the steady simulations, whereas in the unsteady simulations, it defines the time discretization. Both FM and BPF SPL are calculated by taking into account the last five revolutions in the unsteady simulations.
2.2. Aeroacoustics Solver: Farassat Formulation-1A
To capture the tonal noise emitted by MAV rotors, Farassat’s formulation-1A of the Ffowcs Williams and Hawkings acoustic analogy has been implemented in the NVLM code on both unsteady and steady solvers [
20,
79]. The considered acoustic sources are as follows:
The thickness noise depends on the rotor geometry presented in
Figure 3 (through the
normal vector and the integral on the blades surface
), the rotational speed (velocity
and Mach vector
), the position of the observer (
is the position of the observer and
is the observer vector that defines the position of the observer from a radiating point in the geometry), and the properties of the air (density
and speed of sound
c):
with:
The loading noise depends on the loading distribution
l and its derivative
on the blade surface (in the case of the present NVLM code, on the mean camber line of the blade, see
Figure 3). These variables, calculated beforehand by the aerodynamic code, are the input of the following equation:
with
This formulation is applied to solid boundaries (the blades) for the noise calculation. Therefore, only thickness and loading noise are radiated here, while the quadrupole term due to the wake turbulence is neglected.
The NVLM code used in this work calculates the pressure jump between the suction and pressure sides on the blade mean camber line due to the thin airfoil approximation. This pressure difference distribution is then used to compute the loading noise (see
Figure 3).
Since lattices do not radiate to the observer at the same time, the retarded time () must be calculated for each lattice, and the acoustic pressure evaluated at this retarded instant. The contribution of each lattice is then taken into account to calculate the whole rotor acoustic pressure. It is worth noting that for both loading and thickness noises, the near-field terms (characterized by an amplitude decay of ) and far-field terms ( decay) are separated.
Once the temporal signal is calculated, the fast-Fourier-transform (FFT) is computed and the BPF SPL captured. Note that the frequency sampling depends on the previously introduced
step-angle (
Section 2.1). The BPF SPL is the input required by the optimization algorithm.
2.3. Optimization Algorithm and Framework
Optimization techniques play a pivotal role in diverse applications, ranging from engineering to machine learning, finance, and many more [
80]. These techniques focus on finding the best solution(s) within a defined parameter space and by satisfying the constraints set for a given problem (Equation (
11)).
where
is the vector constituted by the design variables that are allowed to change during the optimization. Traditional
gradient-based optimization techniques have proven their effectiveness in searching the optimized parameters in single-objective optimizations.
However, this minimum may be local, rather than global. Gradient-based algorithms are, in fact, very sensitive to the initial guess. If it is located next to a local minimum, this family of algorithms will generally become stuck at the local minimum.
Real-world problems often require multi-objective optimizations, which do not converge to a singular, global optimal solution but rather yield a set of Pareto optimal solutions. This complexity is compounded when the gradient of the objective function is challenging to compute or unavailable. To address this, gradient-free algorithms have been developed, offering a viable approach to calculate Pareto optimal solutions in the absence of gradient information.
In this work, genetic algorithms (GA) are employed to simultaneously optimize two distinct and competing objectives. These algorithms, inspired by Charles Darwin’s theory of natural evolution, are based on the survival of the fittest and on the appearance of crossover combinations and mutations that lead to fitter successive generations. Since they work with a population of solutions, a different set of solutions can be maintained throughout the optimization to reach global minima and/or identify the set of Pareto optimal solutions. The main difference between GAs lies in the survival and selection methods. The NSGA-II [
81] (non-dominated sorting genetic algorithm), already used by the authors in the past [
82], is employed in this work.
To link both aerodynamic and aeroacoustic Python [
83] solvers to the optimization algorithm, the Python-based optimization framework
pymoo version 0.4.2 [
84] is used. It was selected for integrating aerodynamic and aeroacoustic solvers due to its capability for parallel computation, customization of crossover and mutation probabilities and distributions, and the extraction of the Pareto front. The reader is referred to Blank and Deb [
84] for further information.
2.5. Optimization Setup
The optimization described in this paper takes into account multiple manufacturing constraints. The best geometries resulting from the optimization are 3D-printed using stereo-lithography techniques at ISAE-SUPAERO and tested in the anechoic room to verify the performance of the optimized geometries. The starting geometry of the following optimization is the 20 cm diameter rotor previously tested in the ISAE-SUPAERO anechoic room [
9]. The rotor diameter was appositely chosen because it is the maximum diameter that could be 3D-printed in one piece with the available 3D-printer. For this study, only two-bladed rotors were tested because they are easier to balance statically. The NACA0012 aerodynamic coefficients calculated using XFOIL had previously been validated against experimental tests for different Reynolds numbers and angles of attack (as shown in [
76]). This airfoil was therefore appropriately chosen to eliminate a possible source of discrepancy between the tests and the calculations when using lookup tables in the non-linear correction module embedded in the actual NVLM code. Its characteristics and the experimental results (evaluated at 2 N thrust) are summarized in
Table 2.
All of the values in
Table 2 are fixed throughout the optimization except for the chord and pitch distributions. That is, all tested rotors have two blades with a NACA0012 airfoil constant over the entire span of the blade and a diameter of 20 cm.
In order to minimize the number of parameters and still test very different rotor geometries, a blade reconstruction procedure previously introduced by the authors [
82] is used. The optimization algorithm can select the parameter values within predefined minimum and maximum bounds using a continuous function. In this work, pitch and chord distributions are defined by three points: a point located at the root of the blade, one at the tip (called TIP
x, where
x is the parameter), and a control point (CP
x, for which the spanwise position CP
x,pos can be adjusted). The pitch and chord distributions are constructed as follows:
The root pitch angle and chord length are set to 10° and 0.025 m, respectively;
The CPpitch,pos defines the spanwise position at which the value CPpitch,ang is applied. The derivative of the curve at this point is fixed to zero (same strategy for the chord distribution with the CPchord,pos and CPchord,len variables, respectively);
The pitch angle at the tip of the blade is defined by TIPpitch,ang (the chord length at the tip by TIPchord,len).
The
BPoly.from_derivatives function of the Python (version 3.8.2) package
Scipy version 1.4.1 [
86] permits us to obtain the final distribution from the three different points by constructing a piecewise polynomial in the Bernstein basis. The parametrization used in this case has an advantage over conventional NURBS (non-uniform rational basis splines) because the interpolated curve always passes through the control points. This means that both the lower and upper bounds in the GA have a physical meaning.
Figure 5 shows examples of pitch distribution interpolation.
The total design space hence consists of six variables, presented in
Table 3.
As this paper focuses on the performance of optimized 3D-printed rotors, multiple manufacturing and operational constraints are defined and set in the optimization to achieve robust and 3D-printable rotors. These constraints include limits on the chord length and solidity, due to the fact that rotors with higher solidities are generally heavier and more challenging to manufacture. The upper solidity limit was established to match that of the initial rotor design, while the lower limit was set at 0.08, a value that was arbitrarily chosen. Quadcopter drones are easily maneuverable due to their simplicity, as flight control is achieved by changing the rotational speed of the different rotors. However, the higher the inertia of the rotors, the more difficult it is to maneuver the drone. For this reason, an upper limit of inertia is fixed equal to the inertia of the starting geometry. For the rotational speed, because one critical structural mode of the test stand sustaining the rotor in the anechoic room is at around 3000 rpm, this value is selected as the lower limit.
The constraints (summarized in
Table 4) are handled using a technique described by Deb et al. [
81]. In an unconstrained single-objective optimization, when two solutions are compared, the one with the better objective value is chosen. In constrained multi-objective optimization, three different situations can be encountered: both solutions are feasible (i.e., they do not violate any constraints); one is feasible, and the other is not, and both violate the constraints. These three situations lead to different outcomes: the solution with the better objective value is chosen, the feasible solution is chosen, and the solution with smaller overall violation constraint is chosen, respectively. Taking into account solutions not feasible in the generation of a new offspring population is very important during the first generations, when the initial population is randomly chosen and the constraints are tight.
In this study, each generation consists of 100 candidates randomly selected at the beginning. Both crossovers and mutations are activated with a probability of and , respectively, and the optimization stops after 50 generations.