Algorithm for Propeller Optimization Based on Differential Evolution. (2024)

Link/Page Citation

Author(s): Andry Sedelnikov (corresponding author) [*]; Evgenii Kurkin; Jose Gabriel Quijada-Pioquinto; Oleg Lukyanov; Dmitrii Nazarov; Vladislava Chertykovtseva; Ekaterina Kurkina; Van Hung Hoang

1. Introduction

The propeller is one of the most important elements of an aircraft that provides lift and propulsion in the air. The relevance of propellers as a propulsion system has increased with the current growth in the production of unmanned aerial vehicles for various purposes around the world. The air propeller’s design is an important part of aircraft design, beginning from the conceptual design stage [1]. An optimized propeller can significantly reduce emissions, improve acoustic response, and enhance performance [2]. The design of air propellers to minimize propulsion system energy costs began with Zhukovsky, Betz, and Goldstein [3,4] and was then refined by other scientists [5,6,7,8].

Propeller optimization is a complex, often multi-objective problem that requires consideration of many factors and constraints, including aerodynamic issues, strength, and acoustics [7,9,10,11,12]. The most common goals of propeller optimization from an aerodynamic perspective are the requirements of maximum thrust, maximum efficiency, and a balanced propeller based on Paretto optimality [13,14,15].

In recent years, computational fluid dynamics (CFD) using computer-aided engineering (CAE) systems has become an important method for propeller design and computation [16,17,18]. Therefore, the use of the finite volume method is irrational because it requires large computational power. Some of the most commonly used propeller aerodynamic characteristics calculation methods currently are Blade Element Method (BEM), Blade Element Momentum Theory (BEMT), Vortex Lattice Method (VLM), and Computational Fluid Dynamics (CFD) [1,19]. For tip-loss correction, different methods are used [20], including the Goldstein function [21] and Glauert correction [22,23,24]. The BEM and BEMT methods, due to the fast calculation, have been the most used in optimization processes. The Isolated Section Method (ISM) is similar to BEMT method, based on the Zhukovsky propeller vortex theory, using the Goldstein function for tip-loss correction [25,26]. It should be noted that the accuracy of the BEM, BEMT and ISM methods depends on the precision of the propeller’s cross-sections aerodynamic coefficients prediction. XFOIL code [27] based on the discrete vortex method is a common technique for predicting the airfoils aerodynamic characteristics in propeller design [24,28].

The methodology of propeller geometry parametrization is part of its the aerodynamic optimization problem. Bezier surfaces and curves have the possibility of generating complex shapes using some control points [29]. This feature makes the use of Bezier curves and surfaces suitable for describing propeller cross-sections (airfoils [30,31,32]) and the propeller blades in general, including chord and twist [1,33]. Due to the reduced number of control points, these curves and surfaces are used in propeller shape optimization processes [34], especially when used in conjunction with evolutionary algorithms [35,36] and artificial neural networks [1,33].

Various methods are used in solving the propeller optimization problem. The authors of [37] solved a multi-objective optimization problem subject to a set of constraints using a direct search algorithm. The paper [28] describes the experimental verification and application of a method of multi-criteria optimization using genetic algorithms for the design of a propeller for a high-altitude aircraft. The articles [38,39] discuss the effectiveness of differential evolution-based algorithms for optimizing the shape of different types of propellers.

BEMT and IMS methods for calculating a given shape propeller aerodynamic performance require an airfoil angle of attack selection cycle, causing the need for an airfoil drag polar calculation database [40] and drag curve extrapolation when the geometric twist exceeds the limits on angles of attack. Replacing the twist angle in the design variables with the angle of attack allows the optimization process to immediately search for the desired airflow field and eliminate an additional cycle of searching for the angle of attack for each propeller individual. In the current work, a modification of ISM consists in the airfoil effective angle usage instead of the geometrical twist angle is presented. The geometric twist is calculated for the optimum propeller when the optimum distribution of the angle of attack over the span is found. This approach also improves the robustness of the optimization process by directly controlling the constraints imposed on the angle of attack.

The parametric geometric model of the propeller used in the current work implements Bezier curves to specify the dependencies of the parameters of thickness, curvature, and positions of the maximum thickness and maximum curvature of the airfoil along the blade span. This makes it possible to use different thicknesses and curvatures of airfoils along the blade span, which allows the aerodynamic twist of the blade to be varied.

In addition, the article considers the use of interpolation of the calculation of the aerodynamic characteristics of the airfoils on the blade span, which also allows to accelerate the process.

The techniques of parameterization of propeller geometry, calculation of its aerodynamic characteristics and optimization of its shape considered in the paper are implemented in the OpenVINT code [41]. This code uses the following techniques: geometry modeling with Bézier curves, a modified ISM for obtaining propeller aerodynamic coefficients based on XFOIL for airfoils characteristics calculation, optimization by an evolutionary algorithm with adaptive techniques, and population size reduction methods.

The structure of this paper is as follows: Section 2 describes the methods that comprise the OpenVINT algorithm; in Section 3, two case studies are described in which the versatility of the algorithm for optimizing different types of propellers is shown; finally, the discussions relevant to the results and conclusions of this work are presented.

2. Methods

2.1. Mathematical Model of Propeller Optimization

The objective of this paper is to develop an algorithm that allows the optimization of propeller geometry to minimize the power required to achieve a given thrust. This translates into improved energy consumption. Mathematically, the optimization problem can be expressed as a constrained optimization problem. This paper describes the development of an optimization algorithm to improve the performance of propellers used in unmanned aerial vehicles. min W(x), such that Tobj - T(x) = 0, x?X where W(x) is the required power of the propeller; T(x) is the thrust provided by the propeller; T[sub.obj] is the desired thrust; x is the vector of design parameters belonging to the feasible set of solutions X.

For this optimization problem to be solved by evolutionary algorithms, the problem had to be converted into an unconstrained optimization problem. This was achieved by using a penalty function, which will be used as a fitness function. The penalty function is expressed as follows:(1)L(x)={W(x)if?(x)=0R?(x)+U*if?(x)>0?W(x)=U*R?(x)+W(x)if?(x)>0?W(x)>U* where (2)?(x)=max?{0,T[sub.obj]-T(x)} where U* is an upper bound on the constrained global minimum value; R is a penalty parameter that allows to equate the values obtained in ?(x) with U*. A large value is assigned to U* initially, but the value needs to be updated with the current best known function value at feasible points. Hence, the initial U* is not updated until a feasible solution has been found. ?(x) > 0 only if x is infeasible [42]. Therefore, the optimization problem remains as: min L(x), such that x?X

2.2. Selection of Design Variables

Two types of design variables were proposed: those that describe in a general way the geometry of the propeller or the operation, such as the diameter of the propeller (d, in m), number of blades (B), and number of revolutions per minute (n[sub.m], rev/min); and others that change depending on the radius of the propeller, such as the chord length, effective angle of the airfoil, and the geometry of the airfoil. We considered that these variables have a smooth and continuous variation along the propeller radius using Bezier curves. This variation was achieved using two quadratic Bezier curves [43] (see Figure 1). A quadratic Bezier curve is the path drawn by the following function:(3)BC(t)=(1-t[sup.2])P[sub.0]+2(1-t)tP[sub.1]+t[sup.2]P[sub.2] where P[sub.0], P[sub.1] and P[sub.2] are control points, and t is a parameter that always varies from 0 to 1.

BC(t) can also be decomposed as (B[sub.r¯](t), B[sub.x](t)), where r¯ indicates the relative radius of the propeller and x the design parameter. (4)B[sub.r¯](t)=(1-t)[sup.2]r¯[sub.0]+2(1-t)tr¯[sub.1]+t[sup.2]r¯[sub.2], (5)B[sub.x](t)=(1-t)[sup.2]x[sub.0]+2(1-t)tx[sub.1]+t[sup.2]x[sub.2],

The control points that determine the Bézier curves for each design variable as a function of the propeller’s relative radius are as follows:

* Root curve

(6){r¯[sub.0][sup.r]=0.1r¯[sub.1][sup.r]=0.5r¯[sub.xm]+0.05r¯[sub.2][sup.r]=r¯[sub.xm]{x[sub.0][sup.r]=x[sub.r]x[sub.1][sup.r]=x[sub.m]x[sub.2][sup.r]=x[sub.m]

* Tip curve

In Equations (4)–(7), x may be the chord relative to the propeller diameter (c/d), the effective angle (a), maximum thickness (y[sub.t]), position of the maximum thickness relative to the chord (x[sub.t]), the maximum camber (y[sub.c]), and the position of the maximum camber relative to the chord (x[sub.c]). Each of these variables is with respect to the profile of each section of the propeller blade. The variable r¯[sub.xm] refers to the union position of the two Bezier curves. The subscripts r, m, t refers to the values of the design variable with respect to 10% of the blade radius, the union point of the Bezier curves, and 97% of the blade radius, respectively. Therefore, each vector of design variables is formed as x[sub.i] = [(c/d)[sub.r], (c/d)[sub.m], (c/d)[sub.t], r¯[sub.cm], a[sub.r], a[sub.m], a[sub.t], r¯[sub.am], y[sub.tr], y[sub.tm], y[sub.tt], r¯[sub.ytm], x[sub.tr], x[sub.tm], x[sub.tt], r¯[sub.xtm], y[sub.cr], y[sub.cm], y[sub.ct], r¯[sub.ycm], x[sub.cr], x[sub.cm], x[sub.ct], r¯[sub.xcm], n[sub.m], B, d][sub.i].

2.3. Development of the Input Geometry

The Isolated Sections Method (ISM) [44,45] was selected to obtain the aerodynamic characteristics of the different propeller configurations. This method requires knowledge of the geometry, chord, and effective angle of attack of the airfoil in different sections of the propeller blade. The chord values and the effective angle of attack were obtained by constructing the distribution curves (Bezier curve) and using the input parameters (c/d)[sub.r], (c/d)[sub.m], (c/d)[sub.t], r¯[sub.cm], a[sub.r], a[sub.m], a[sub.t], and r¯[sub.am]. To obtain the coordinates of the airfoil in each section, an airfoil parameterization method based on the Bezier-PARSEC technique was proposed [46]. The use of Bezier curves makes it possible to model a wide variety of profile types [39,47]. Each airfoil is constructed with four cubic Bezier curves, two curves for determining the thickness of the airfoil and two curves for determining the camber of the airfoil (see Figure 2).

The parametric functions for determining each cubic Bezier curve are as follows:(8)B[sub.x](t)=(1-t)[sup.3]x[sub.0]+3(1-t)[sup.2]tx[sub.1]+3(1-t)t[sup.2]x[sub.2]+t[sup.3]x[sub.3], (9)B[sub.y](t)=(1-t)[sup.3]y[sub.0]+3(1-t)[sup.2]ty[sub.1]+3(1-t)t[sup.2]y[sub.2]+t[sup.3]y[sub.3],

The construction of the profile is carried out using the following equations:(10)X[sub.t]=B[sub.xt][sup.le](t)+B[sub.xt][sup.te](t) (11)Y[sub.t]=B[sub.yt][sup.le](t)+B[sub.yt][sup.te](t) (12)X[sub.c]=B[sub.xc][sup.le](t)+B[sub.xc][sup.te](t) (13)Y[sub.c]=B[sub.yc][sup.le](t)+B[sub.yc][sup.te](t) (14)?=tan[sup.-1]?(dYc/dXc)

The control points for the leading edge thickness curve are defined by:(15){x[sub.0][sup.le]=0x[sub.1][sup.le]=0x[sub.2][sup.le]=0.5x[sub.t]x[sub.3][sup.le]=x[sub.t]{y[sub.0][sup.le]=0y[sub.1][sup.le]=0.34y[sub.t]y[sub.2][sup.le]=0.5y[sub.t]y[sub.3][sup.le]=0.5y[sub.t]

The control points for the trailing edge thickness curve are defined by:(16){x[sub.0][sup.te]=x[sub.t]x[sub.1][sup.te]=0.3+0.7x[sub.t]x[sub.2][sup.te]=0.6+0.4x[sub.t]x[sub.3][sup.te]=1{y[sub.0][sup.te]=0.5y[sub.t]y[sub.1][sup.te]=0.5y[sub.t]y[sub.2][sup.te]=0.29y[sub.t]y[sub.3][sup.te]=0

The control points for the leading edge camber curve are defined by:(17){x[sub.0][sup.le]=0x[sub.1][sup.le]=x[sub.c]/3x[sub.2][sup.le]=2x[sub.c]/3x[sub.3][sup.le]=x[sub.c]{y[sub.0][sup.le]=0y[sub.1][sup.le]=0.71y[sub.c]y[sub.2][sup.le]=y[sub.c]y[sub.3][sup.le]=y[sub.c]

And the control points for the trailing edge camber curve are defined by:(18){x[sub.0][sup.te]=x[sub.c]x[sub.1][sup.te]=(1+2x[sub.c])/3x[sub.2][sup.te]=(2+x[sub.c])/3x[sub.3][sup.te]=1{y[sub.0][sup.te]=y[sub.c]y[sub.1][sup.te]=y[sub.c]y[sub.2][sup.te]=0.43y[sub.c]y[sub.3][sup.te]=0

The points for determining the upper curve of the airfoil are determined by:(19)X[sub.U]=X[sub.c]-Y[sub.t]sin?? (20)Y[sub.U]=Y[sub.c]+Y[sub.t]cos??

And the points for determining the lower curve of the airfoil are determined by:(21)X[sub.L]=X[sub.c]+Y[sub.t]sin?? (22)Y[sub.L]=Y[sub.c]-Y[sub.t]cos??

Algorithm 1 shows the procedure for obtaining the input conditions for the ISM.

Algorithm 1: Subroutine for creation of the Bezier curves and airfoils

Inputs: x[sub.i], d, r¯, NS

Outputs: (c/d)[sub.i], a[sub.i], x[sub.ti], y[sub.ti], x[sub.ci], y[sub.ci], X[sub.U], Y[sub.U], X[sub.L], Y[sub.L]

1

Create the distribution curve for (c/d)[sub.i] as a function of the r¯[sub.i] of the blade with (4), (5), (6), (7);

2

Create the distribution curve for a[sub.i] as a function of the r¯[sub.i] of the blade with (4), (5), (6), (7);

3

Create the distribution curve for x[sub.ti] as a function of the r¯[sub.i] of the blade with (4), (5), (6), (7);

4

Create the distribution curve for y[sub.ti] as a function of the r¯[sub.i] of the blade with (4), (5), (6), (7);

5

Create the distribution curve for x[sub.ci] as a function of the r¯[sub.i] of the blade with (4), (5), (6), (7);

6

Create the distribution curve for y[sub.ci] as a function of the r¯[sub.i] of the blade with (4), (5), (6), (7);

//Get the airfoil in each section of the blade

7

for s = 1 to NS do

8

Get x[sub.t](r¯[sub.s,i]), y[sub.t](r¯[sub.s,i]), x[sub.c](r¯[sub.s,i]) and y[sub.c](r¯[sub.s,i]);

9

Get X[sub.ts,i] and Y[sub.ts,i] with (8), (9), (10), (11), (15), (16);

10

Get X[sub.cs,i] and Y[sub.cs,i] with (8), (9), (12), (13), (17), (18);

11

Get ?[sub.s,i] with (14);

12

Get the points of the sth-airfoil X[sub.U], Y[sub.U], X[sub.L], Y[sub.L] with (19), (20), (21), (22);

13

return Outputs

2.4. Construction of the Output Geometry

The final geometry of the blade is obtained after executing the ISM algorithm because this process is where the blade twist is obtained. Figure 3 shows the process of drawing a propeller blade from the proposed Bezier curves and the geometric torsion distribution obtained from Algorithm 2. To create the base shape of the blade, the torsion blade axis is first defined, which crosses the chord of each station at the point where the maximum thickness of the aerodynamic profile (x[sub.t]) is located. From this axis, the strings are rotated according to the curve f. On the other hand, with the four Bezier curves x[sub.t], y[sub.t], x[sub.c], and y[sub.c], the aerodynamic profiles of each station are generated using the procedure shown in Figure 2. Finally, the profiles are located and scaled according to the corresponding chord length.

2.5. Isolate Section Method

The objective function of the optimization process is to calculate the propeller power (W, in W) required to achieve thrust (T, in N). A modified version of was used to obtain these values. The variation between the original and our modified version of the isolated section method (ISM) is the geometric twisting of the sections (f, in deg) using the effective angle of attack of the sections (a, in deg) as the input value (Algorithm 2).

The input variables required by the method are as follows: B, V[sub.8], d, n (n, rev/s), and the physical characteristics of the air (density (?, in kg/m[sup.3]), kinematic viscosity (?, [m[sup.2]/s]), and sound speed (a, [m/s])). The isolated section method requires sectioning one of the propeller blades into a finite number of sections (NS). In each section, it is necessary to know the following values: chord length (c, in m), r¯, and the geometry of the airfoil ((X[sub.U], Y[sub.U]), (X[sub.L], Y[sub.L]), where each coordinate is in m).

The iterative method for obtaining the thrust provided by the propeller and the required power is described in [44]. This method was modified to use the effective angle of attack (a) of the profile in each section as a constant and the geometric twist (f) of the section as a variable that is updated in each iteration.

Having the local characteristics of the flow, in addition to the chord, geometry, and angle of attack of the airfoil, coefficients of lift (c[sub.l]) and drag (c[sub.d]) coefficients in each section are calculated. The next step is to obtain U¯[sub.1], V¯[sub.1] and G¯[sub.r] by an iterative process, which is necessary to obtain the coefficients c[sub.t] and m[sub.k]. The equations involved in this process are as follows:(23)v¯[sub.1]=-v¯/2+[square root of v¯2/4+u¯[sub.1](r¯-u¯[sub.1])+2?r¯1u¯12/r¯dr¯], (24)U¯[sub.1]=r¯-u¯[sub.1], (25)V¯[sub.1]=v¯+v¯[sub.1], (26)W¯[sub.1]=[square root of V¯[sub.1][sup.2]+U¯[sub.1][sup.2]], (27)ß[sub.1]=tan[sup.-1]?(V¯1/U¯1), (28)s=Bc/Rp, (29)G¯[sub.r]=1/8sc[sub.l]W¯[sub.1], (30)f[sub.r]=2/pcos[sup.-1]?(e[sup.-0.5B(1-r¯)r¯sin?ß1]).

Previously before initializing the iterative process, it is important to initialize the n values of u¯[sub.1] and ?[sub.r¯][sup.1]u¯12/r¯dr¯ equal to zero. At the end of each iteration u¯[sub.1] has to be updated making use by:(31)u¯[sub.1]=G¯r/frr¯, while the values of ?[sub.r¯][sup.1]u¯12/r¯dr¯ are updated as shown in lines 24, 25 and 26 of Algorithm 2, where trapz (Y, X) integrates along the given axis using the compound trapezoidal rule, Y is the input array to integrate and X is the sample points corresponding to the Y values.

It has been proven that the loop described between lines 10 and 26 of Algorithm 2 only needs 10 cycles to obtain good results.

To obtain the coefficients c[sub.t] and m[sub.k], it is first necessary to calculate their differentials in each section, and then integrate with respect to r¯, making use of trapz(). (32)dc[sub.t]=8G¯[sub.r](U¯[sub.1]-K[sup.-1]V¯[sub.1]) (33)dm[sub.k]=8G¯[sub.r](V¯[sub.1]+K[sup.-1]U¯[sub.1])r¯

By obtaining the coefficients c[sub.t] and m[sub.k], the thrust provided by the propeller (T) and the required power of the propeller (W) can be calculated. (34)T=0.5c[sub.t]??R[sup.2]pR[sup.2] (35)W=0.5m[sub.k]??R[sup.3]pR[sup.2]

Other propeller performance metrics that can be obtained with this method are the dynamic efficiency ?[sub.d] and the static efficiency ?[sub.s] of the propeller. (36)a[sub.p]=T/?n2d4 (37)ß[sub.p]=W/?n3d5 (38)?[sub.p]=V8/nd (39)?[sub.d]=ap?p/ßp (40)?[sub.s]=ct3/2/2mk

Algorithm 2: Modified Isolate Section Method

Inputs: n, B, V[sub.8], d, ?, ?,a,c,r¯, [X[sub.U], Y[sub.U], X[sub.L], Y[sub.L]], NS

Outputs: T, W, ?[sub.d], ?[sub.s]

1

R = d/2;

2

Get ? with (31);

3

Getv¯ with (32);

4

for s = 1 to NS do

5

r[sub.s]=r¯[sub.s]R;

6

Get Re[sub.s] and M[sub.s] in each section;

7

c[sub.ls], c[sub.ds] = runXFOIL(a[sub.s], [X[sub.U], Y[sub.U], X[sub.L], Y[sub.L]][sub.s], Re[sub.s], M[sub.s]);

8

0?u¯[sub.1s],0??r¯1u¯12r¯dr¯[sub.s]

9

for t = 1 to 10 do

10

I[sub.u] = Ø;

11

for s = 1 to NS do

12

Getv¯[sub.1s] with (23);

13

GetU¯[sub.1s] with (24);

14

GetV¯[sub.1s] with (25);

15

GetW¯[sub.1s] with (26);

16

Get ß[sub.1s] with (27);

17

f[sub.s] = a[sub.s] + ß[sub.1s];

18

K[sub.s] = c[sub.ls]/c[sub.ds];

19

Get s[sub.s] with (28);

20

GetG¯[sub.rs] with (29);

21

Get f[sub.rs] with (30);

22

Updateu¯[sub.1s] with (31);

23

u¯1s2r¯s[sub.s]?I[sub.us];

24

for s = 1 to NS do

25

Get?r¯1u¯12r¯dr¯[sub.s]withtrapz(I[sub.u][s:end],r¯[s:end]);

26

for s = 1 to NS do

27

Getdc[sub.ts] with (32);

28

Getdm[sub.ks] with (33);

29

c[sub.t]=trapz(dc[sub.t],r¯);

30

m[sub.k]=trapz(dm[sub.k],r¯);

31

Get a[sub.p], ß[sub.p] and ?[sub.p] with (36), (37) and (38) respectively;

32

Get outputs T, W, ?[sub.d] and ?[sub.s] with (34), (35), (39) and (40) respectively;

33

Save f;

34

return Outputs

2.6. Airfoil Aerodynamic Coefficients Calculation

The aerodynamic coefficients of each section of the blade was obtained by using the Xfoil program, which shows a good performance for the evaluation of profiles. In Figure 3, the experimental coefficients of the Clark-Y profile can be observed using XFOIL, OpenFOAM (with the k-? SST turbulence model) and experimental tests.

In Algorithm 2, the modified ISM is described in pseudo code. To simplify the ISM calculations, it is necessary to obtain the radius (R, [m]), the angular velocity (?, [rad/s]) and the relative velocity (v¯, dimensionless) of the propeller. (41)?=2pn (42)v¯=V8/?R

Sections are defined by the r¯ array. For example, in the cases that we evaluate, the sections begin to be numbered from 10% of the R of the propeller and end at 97% of R (r¯ = [0.1, …, 0.97]). In our modified ISM, having the a of each section as input, it is necessary to calculate the aerodynamic coefficients of each section delimited by r¯. For this, it is necessary to obtain the local Reynolds (Re) and Mach (M) numbers.

To determine the aerodynamic characteristics of the propeller and reduce the number of calculations, interpolation of the profile data was performed over a smaller number of sections, in which the characteristics of the angle of attack were calculated considering the local Reynolds and Mach numbers (see Figure 4).

A study was conducted out on the required number of sections, which showed that the use of 15 sections to calculate the aerodynamics of the profile and 75 sections to integrate thrust and moment allows one to quickly obtain accurate values of the propeller characteristics (Figure 5).

2.7. Isolate Section Method Validation

The proposed method of propeller calculation was carried out by comparing it with the experimental data of the NACA 5868–9 propeller (experimental results, geometric and kinematic characteristics of the propeller are given in [45]), as well as by comparison with the results of numerical mathematical modeling by solving the Navier–Stokes equations in ANSYS CFX 18.2 software using the shear stress transport turbulence model. The propeller diameter is 3.15 m. The problem statement and the dimensions of the static and rotational domains are shown in Figure 6. The configuration and dimensions of the domain are based on those shown in [48]. The stator mesh has 2,618,015 tetrahedral elements, and the rotor mesh has 11,130,389 elements (tetrahedra and pyramids). The minimum element sizes were 150 mm for the stationary domain, 15 mm for the rotating domain, 4.0 mm on the propeller wall. The first boundary layer thickness was 0.05 mm. In the CFX calculation, the propeller was modeled together with the nacelle, matching the geometry to the experimental conditions. The free wind speed is 170 km/h, with rotational speeds from 800 rpm to 2500 rpm. Validation was performed for two cases of propeller geometry, with blade angles of f[sub.0.75] = 25° and f[sub.0.75] = 35°. The values of the mean residuals for all equations do not exceed 10[sup.-3]. The velocity field and pressure distribution are shown in Figure 7. The comparison results in terms of thrust coefficient and power from the advanced ratio are shown in Figure 8.

2.8. Optimization Algorithm

The optimization algorithm used to solve the task is based on the Differential Evolution (DE) algorithm. The proposed algorithm contemplates the use of self-adaptive schemes of evolutionary operators, methods of Population Size Reduction (PSR), sampling techniques for the selection of individuals from the initial population, and stopping conditions that adapt to the needs of the optimization process. In addition, the developed algorithm contemplates the use of parallel computing strategies to accelerate the calculation of the values of the objective function.

Differential evolution is a stochastic, population-based algorithm developed for real-valued function optimization problems. It operates by having a population of individuals (x vector) that move around in the search space by recombining through crossover and mutation with other existing individuals in the population. Through a selection process, a newly generated individual is accepted as part of the population if the new individual x is an improvement; otherwise, it is discarded. This iteration process is repeated to find a vector x that optimizes a function f(x) [49]. As with other evolutionary algorithms, the search performance of differential evolution algorithms depends on the control parameter settings. A standard differential evolution has three main control parameters: population size, scaling factor F, and crossover factor CR. However, it is well-known that the optimal settings of these parameters are problem dependent. Therefore, when applying differential evolution to a real-world problem, it is often necessary to tune the control parameters to obtain the desired results. In practical cases, many researchers suggest the use of self-adaptive schemes to adjust online control parameters during the search process. A variant of differential evolution that applies this type of scheme is the success history-based adaptation for differential evolution (SHADE) [50]. In previous works by the authors of this article, the good performance of the SHADE algorithm for the optimization of aerodynamic bodies such as wings has been corroborated [51].

In SHADE, the mutation factor F ? [0, 1] controls the magnitude of the differential mutation operator and CR ? [0, 1] is the crossover factor. In the beginning, the contents of MCR, k, MF, k (k = 1, …, H) are all initialized to 0.5. In each generation g, the control parameters CRi and Fi used by each individual xi are generated by randomly selecting an index ri from [1, H], and then applying the formulas below:(43)CR[sub.i]={0ifM[sub.CR,ri]=-1randn[sub.i](M[sub.CR,ri],0.1)otherwise (44)F[sub.i]=randc[sub.i](M[sub.F,ri],0.1)

Here, randn[sub.i](M, 0.1), randc[sub.i](M, 0.1) are values selected randomly from normal and Cauchy distributions with mean M and variance 0.1. In case a value for CR[sub.i] outside of [0, 1] is generated, it is replaced by the limit value (0 or 1) closest to the generated value. When F[sub.i] > 1, F[sub.i] is truncated to 1, and when F[sub.i] = 0, Equation (44) is repeatedly applied to try to generate a valid value. In Equation (43), if M[sub.CR,ri] has been assigned the “terminal value” -1, CR[sub.i] is set to 0.

The mutation strategy used by SHADE current-to-pbest/1 is a variant of the current-to-best/1 strategy where the greediness is adjustable using a parameter p. (45)v[sub.i,g]=x[sub.i,g]+F[sub.i](x[sub.pbest,g]-x[sub.i,g])+F[sub.i](x[sub.r1,g]-x[sub.r2,g])

The individual x[sub.pbest,g] is randomly selected from the top p·NP members in the g-th generation (p ? [0, 1]). F[sub.i] is the F parameter used by individual x[sub.i]. The greediness of current-to-pbest/1 depends on the control parameter p to balance exploitation and exploration (small p behaves more greedily).

SHADE uses binary crossover strategy, (46)u[sub.j,i,g]={v[sub.j,i,g]ifrand(0,1)=CR[sub.i]x[sub.j,i,g]otherwise

Perform selection between the trial vector (u[sub.i,g]) and target vector (x[sub.i,g]) using a greedy selection criterion, (47)x[sub.i,g+1]={u[sub.i,g]iff(u[sub.i,g])=f(x[sub.i,g])x[sub.i,g]otherwise

SHADE uses a historical memory M[sub.CR], M[sub.F] which stores a set of CR, F values that have performed well in the past, and generates new CR, F pairs by directly sampling the parameter space close to one of these stored pairs.

In Algorithm 3, index k determines the position in memory to be updated. In generation g, the k-element in the memory is updated. At the beginning of the search, k is initialized to 1. k is incremented whenever a new element is inserted into the history. If k > H, k is set to 1. In the updated algorithm 1, note that when all individuals in generation g fail to generate a trial vector that is better than the parent, i.e., S[sub.CR] = S[sub.F] = Ø, the memory is not updated. The weighted Lehmer mean mean[sub.WL](S) is computed using the following formula:(48)mean[sub.WL](S)=?m=1SwmSm2/?m=1SwmSm (49)w[sub.m]=?fm/?l=1S?fl (50)?f[sub.m]=|f(u[sub.m,g])-f(x[sub.m,g])|

Algorithm 3: Memory update algorithm in SHADE

Inputs: S[sub.CR], S[sub.F], M[sub.CR,k,g], M[sub.F,k,g], k, H

Outputs: M[sub.CR,k,g+1], M[sub.F,k,g+1], k

1

if S[sub.CR] ? Ø and S[sub.F] ? Ø then

2

if M[sub.CR,k,g] = -1 or max(S[sub.CR]) = 0 then

3

M[sub.CR,k,g+1] = -1;

4

else

5

M[sub.CR,k,g+1] = mean[sub.WL](S[sub.CR]);

6

M[sub.F,k,g+1] = mean[sub.WL](S[sub.F]);

7

k++;

8

if k > H then

9

k = 1;

10

else

11

M[sub.CR,k,g+1] = M[sub.CR,k,g];

12

M[sub.F,k,g+1] = M[sub.F,k,g];

13

return Outputs

The amount of fitness improvement ?f[sub.m] is used to influence the parameter adaptation (S refers to either S[sub.CR] or S[sub.F]). As M[sub.CR] is updated, if M[sub.CR,k,g] = -1 or max(S[sub.CR]) = 0 (i.e., all elements of S[sub.CR] are 0), M[sub.CR,k,g+1] is set to -1. Thus, if M[sub.CR] is assigned the terminal value -1, then M[sub.CR] will remain fixed at -1 until the end of the search. This has the effect of locking CR[sub.i] to 0 until the end of the search, causing the algorithm to enforce a “change-one-parameter-at-a-time” policy, which tends to slow down convergence, and is effective on multimodal problems.

The SHADE algorithm works well in conjunction with the PSR methods. To develop this optimization algorithm, we decided to incorporate the continuous adaptive population reduction (CAPR) method. The CAPR method gradually reduces the population size according to the change in the gradient of the fitness value [52]. (51)NP[sub.g+1]={?[sub.g]/?[sub.g-1]?0<?[sub.g]/?[sub.g-1]<1NP[sub.g]otherwise (52)NP[sub.g+1]={NP[sub.g+1]NP[sub.g+1]>NP[sub.min]NP[sub.min]NP[sub.g+1]=NP[sub.min] where (53)?[sub.g]=favgxg-favgxg-1/favgxg,?[sub.g-1]=favgxg-1-favgxg-2/favgxg-1

In the third generation and in subsequent generations, the evaluated function values of all vectors in the population are averaged to be f[sub.avg](x[sub.g]). This value, together with that form the previous evaluation generation, is used to calculate the normalized gradient value ?[sub.g]. ?[sub.g-1] is calculated in a similar fashion using the previous average function evaluation value, f[sub.avg](x[sub.g-1]), and the one before the previous f[sub.avg](x[sub.g-2]). If the ratio ?[sub.g]/?[sub.g-1] is within the range of [0, 1], then NP is reduced by a fraction equal to the ?-th root of the ratio ?[sub.g]/?[sub.g-1]. The reason for taking root of the ratio is to slow down the population size reduction rate.

Another criterion to consider when increasing the performance of algorithms based on differential evolution is the generation of the initial population. It has been shown that a population, whose individuals are best distributed throughout the entire design space, has a greater chance of finding a global optimum, in addition to reducing the search time. The LHS design is a statistical method for generating a quasi-random sampling distribution. It is one of the most popular sampling techniques used in computer experiments because of its simplicity and projection properties with high-dimensional problems. The LHS is built as follows: each dimensional space, representing a variable, is cut into n sections, where n is the number of sampling points, and only one point is placed in each section [53].

For real-case optimization processes, it is common to use two types of stopping criteria. The following two types of stopping criteria were considered for this algorithm [54]:

* Exhaustion-based criteria: due to limited computational resources, the optimization run might be terminated after a certain number of objective function evaluations or CPU time. Typically, a maximum number of generations or number of objective function evaluations is used in combination with every stopping criterion to prevent the algorithm from running forever if a criterion cannot stop the run.

* Distribution-based criteria: for differential evolution algorithms, all individuals eventually converge to the optimum. Therefore, it can be concluded that convergence occurs when individuals are close to each other. Because it is assumed that the optimum is not known for the reference criterion, the distances between the population members are examined. This type of criterion can be applied to the design space or objective space.

One of the main disadvantages of evolutionary algorithms is that they need to evaluate multiple vectors to find the global optimum, which implies calculating the values of the objective function many times. One way to speed up the calculation process is to use parallel computing strategies. The strategy that was proposed to be used is Lightweight Pipelining (LP) [55]. The pipelining process provides an easy approach for downloading and using the models on demand. It helps in parallelization, which means that different jobs can be run in parallel. It also reduces redundancy and helps to inspect and debug the data flow in the model. Some features that pipeline provides are on-demand computing, tracking of data and computation, and inspection of data flow.

OpenVINT is the union and adaptation of each of the aforementioned algorithms and methods to achieve the objective of the optimization process mentioned at the beginning of this work. Coding of this algorithm was performed mainly on Python 3 in a GNU/LINUX environment.

3. Results

3.1. Case Studies

To evaluate the performance of the OpenVINT algorithm, two cases were performed:

* case 1—obtain the optimal design of a propeller used for an engine of a fixed-wing aircraft, considering a free flow velocity of 25 m/s and a minimum required thrust is 7.5 N;

* case 2—obtain the optimal design of a propeller used in a helicopter, considering a free flow velocity of 2 m/s and a minimum thrust required is 6.5 N.

In both cases, the following physical properties of the air flow were used: ?, 1.225 kg/m[sup.3]; ?, 0.000014607 m[sup.2]/s; and a, 340.294 m/s. The intervals of the design variables implemented in each of the tests are shown in Table 1.

For both optimization cases the following parameters of the optimization algorithm were taken:

* in real optimization problems it is considered to use at least 50 individuals in the initial population [56], and as a minimum population 10 individuals were considered;

* the stop conditions contemplated were a maximum number of evaluated generations of 200, and an e value of 1 W to fulfill the condition indicated in line 45 of Algorithm 4;

* a ? factor of 50 was used in (51);

* finally, an initial U* value of 350 W was considered.

All the tests were performed on a PC whose characteristics are as follows:

* operative system, Ubuntu 22.04.3 LTS;

* processor, 13th Gen Intel[sup.®] Core[sup.TM] i5-13400Fx12;

* RAM memory, 64 GB;

Algorithm 4: OpenVINT algorithm

Inputs: d, V[sub.8], ?, ?, a, T[sub.min], r¯, NS, [X[sub.min], X[sub.max]], G, NP, NP[sub.min], e, U*, ?, H, p

Outputs: x[sub.opt], L(x[sub.opt])

//Initialization phase

1

g = 1;

2

Initialize of metrics;

3

Initialize population P[sub.g] with LHS;

//Parallelized loop by joblib

4

for i = 1 to NP do

5

Apply Algorithm 1 for x[sub.i,g];

//Parallelized loop by joblib

6

for i = 1 to NP do

7

Get T(x[sub.i,g]), W(x[sub.i,g]), ?[sub.d](x[sub.i,g]), ?[sub.s](x[sub.i,g]) with Algorithm 2;

8

Get ?(x[sub.g]) with (2) and L(x[sub.g]) with (1);

9

Update U*;

10

Save L[sub.avg](x[sub.g]);

11

Save data of generation g;

12

Set all values in M[sub.CR], M[sub.F] to 0.5;

13

Archive A = Ø;

14

k = 1;

//Main loop

15

for g = 2 to G do

16

S[sub.CR] = Ø, S[sub.F] = Ø;

17

for i = 1 to NP do

18

r[sub.i] = select from [1, H] randomly;

19

Get CR[sub.i,g] with (43);

20

Get F[sub.i,g] with (44);

21

Get mutation vector v[sub.i,g] with (45);

22

Get trial vector u[sub.i,g] with (46);

//Parallelized loop by joblib

23

for i = 1 to NP do

24

Apply Algorithm 1 with u[sub.i,g];

//Parallelized loop by joblib

25

for i = 1 to NP do

26

Get T(u[sub.i,g]), W(u[sub.i,g]), ?[sub.d](u[sub.i,g]), ?[sub.s](u[sub.i,g]) with Algorithm 2;

27

Get ?(u[sub.g]) with (2) and L(u[sub.g]) with (1);

28

Update U*;

29

for i = 1 to NP do

30

if L(u[sub.i,g]) = L(x[sub.i,g]) then

31

x[sub.i,g+1] = u[sub.i,g];

32

x[sub.i,g]?A;

33

CR[sub.i,g]?S[sub.CR], F[sub.i,g]?S[sub.F];

34

else

35

x[sub.i,g+1] = x[sub.i,g];

36

Update memories M[sub.CR] and M[sub.F] with Algorithm 3;

37

Save L[sub.avg](x[sub.g+1]);

38

if g = 3 then

39

Get ?[sub.g] and ?[sub.g-1] with (53);

40

Get NP[sub.g+1] with (51);

41

if NP[sub.g+1] < NP[sub.min]then

42

Apply (52);

43

(NP[sub.g] - NP[sub.g+1])-th worst vectors ? A;

44

Save data of generation g+1;

45

if |L[sub.avg](x[sub.g+1]) - L[sub.opt](x[sub.g+1])| =e then

46

break;

47

k++;

48

Print metrics plots;

49

Output

50

Drawing the optimal propeller in point clouds;

3.2. Optimization Results

For case 1, the following parameters of the optimization algorithm were taken:

* 6 tests were performed, 3 of them were with an initial population of 50 vectors and the other 3 with an initial population of 100 vectors;

* a minimum population of 10 vectors was considered in all tests;

* the stop conditions contemplated were a maximum number of evaluated generations of 200, and an e value of 1 W to fulfill the condition indicated in line 48 of algorithm 4;

* a ? factor of 50 was used in (51);

* finally, an initial U* value of 350 W was considered.

Figure 9 shows the convergence graphs of the metrics used to evaluate the propeller optimization process for case 1.

Table 2 shows the values of the optimal design vector as well as the optimal values of the required power output, thrust provided, and the dynamic efficiency of the propeller. In addition, the average values, and the coefficient of variation (CV) for each parameter are shown, which allows evaluation of the hom*ogeneity of the data.

The geometric differences between each of the blades obtained for case 1 can be visualized in Figure 10 and Figure 11.

For case 2 the following parameters of the optimization algorithm were taken:

* A total of 3 tests were performed, with an initial population of 50 vectors and a minimum population of 10 vectors;

* the stop conditions contemplated were a maximum number of evaluated generations of 200, and an e value of 1 W to fulfill the condition indicated in line 48 of algorithm 4;

* a ? factor of 50 was used in (51);

* finally, an initial U* value of 350 W was considered.

Figure 12 shows the convergence graphs of the metrics used to evaluate the propeller optimization process. The final values of the optimization are shown in more detail in Table 3.

Graphically, the differences between the blades obtained in each of the tests can be observed in Figure 13 and Figure 14.

3.3. Experimental Validation

To verify the results obtained in the optimization process, one of the propellers obtained in the tests performed for case 2 (the propeller obtained in test 3) was selected. The propeller was manufactured by contact molding continuous carbon fibers using epoxy thermosetting binder in a closed mold. The mold was fabricated using a CNC milling machine from the geometry obtained in test 3 (Figure 15). The outer layers and the spars passing through the sleeve comprise unidirectional fibers, and the inner two layers are stacked from plain fabric with fiber orientation ±45°.

The proposed experimental setup allows the monitoring of the angular velocity, thrust, and torque of the propeller. The angular velocity of the propeller was monitored using two different methods: measurement from the electric motor with a PWM frequency meter and measurement from the propeller with a laser photo tachometer. The thrust and torque of the propeller are measured using strain gauges at the base of the motor. The power requirement of the propeller was calculated as the product of torque and angular velocity (see Figure 16) [57].

A comparison was made of the dependences of thrust on the required power for two propellers—designed (test 2, case 3) and commercial DJI 10 × 4.5'' propeller (Figure 17). Seven experimental tests (T1…T7) were carried out for each propeller, varying in rotational speed (Figure 18). Each test consists of from 6 to 8 measurements—a total of 46 measurements for designed propeller and 49 measurements for DJI propeller. The dependences of thrust on required power for each test were interpolated by smooth univariable spline of 3 degrees and then lines of average values and scatter equal to the standard deviation are plotted (Figure 19).

Comparison of calculated by ISM and experimental data for the designed propeller (Figure 19a) showed that for the main propeller regimes (trust more than 2N) the difference does not exceed 3.5%. A comparison between the designed propeller and a commercially available propeller (DJI 10 × 4.5'' propeller) is shown in Figure 19b. The proposed algorithm demonstrates that it can design a propeller that provides a higher thrust (from 10 to 15%) with the same power requirement as commercial propellers in main operating mode.

4. Discussion

In the results of the tests carried out, presented in Table 2 and Table 3, the following points could be observed:

* for the optimization of a helicopter propeller, it is sufficient to use an initial population of 50 vectors to obtain similar results (CV values less than 30%);

* to optimize a tractor propeller, it was necessary to use an initial population of 100 vectors to improve the search for a global optimal;

* in both cases, there is a greater variation in the parameters that describe the geometry of the blade root, whereas in the parameters that describe the rest of the blade geometry, the results were more hom*ogeneous.

Figure 9 and Figure 12 illustrate that in the most representative blade section (in 75% of the blade length), practically similar cross-sectional (airfoil) geometries are obtained. Similar aerodynamic behavior implies that similar values of thrust provided and thrust required are obtained.

When performing the experimental tests with one of the propellers obtained in case 2, it was observed that the aerodynamic characteristics of the propeller obtained with the help of the ISM and XFOIL were comparable to those obtained in the experimental tests, as shown in Figure 19a. These results provide certainty in the calculated values of the objective function for each vector evaluated in the optimization process. While in the comparison that was made with one of the optimized propellers and with a commercial propeller, it was observed that the OpenVINT algorithm provides new propeller options that require less power to achieve a certain thrust.

Regarding the calculation speed of the algorithm, it can be observed that the OpenVINT algorithm provides better performance when optimizing propellers used in helicopters. For this type of propeller, the algorithm requires approximately 1 h of calculation (using a computer as described in Section 3.1), whereas optimizing a tractor propeller for the engine of a fixed-wing aircraft requires approximately 3 h of calculation. Of course, these times can still improve the design space for each type of propeller. This type of analysis is proposed for future studies.

5. Conclusions

The paper presents an algorithm for propeller shape optimization based on isolated section method, similar to BEMT. Using Bezier curves to specify the dependencies of the parameters of thickness, curvature, and positions of the maximum thickness and maximum curvature of the airfoil along the blade span allows the aerodynamic twist of the blade to be varied.

The usage of a cross-section airfoil angle of attack as a function of blade span as design variable adds robustness to the optimization algorithm, protecting from calculation airfoil characteristics at supercritical angles of attack and reduce the number of calculations during the optimization. This makes it possible to reduce the optimization time and increase the accuracy of the optimization solution with geometric models with a large number of design variables and cross-sections used for calculations. Interpolation of the aerodynamic characteristics of sections along the blade span also makes it possible to reduce calculation time. The proposed algorithm allows for obtaining the optimal geometry of the pusher and tractor propellers based on 27 design variables, containing a variable along the blade span airfoil, a nonlinear distribution of twist and chord along the span in three hours of calculation on a personal computer.

In the tests to which OpenVINT optimization algorithm was subjected, it was observed that the parameterization of the proposed propeller geometry is good but still has room for improvement. From the 27 proposed design parameters, the parameters describing the shape of the blade root have little importance in determining the aerodynamic performance of the propeller, which implies that the number of design parameters can be reduced, or the value of these parameters can be refined using multidisciplinary optimization, for example based on stiffness and strength constraints.

The experimental tests validate the calculation method and confirm that the algorithm in general can provide configurations with better aerodynamic performance to the propellers that are available on the market.

Author Contributions

Conceptualization, D.N., E.K. (Evgenii Kurkin) and O.L.; methodology, D.N., J.G.Q.-P., E.K. (Evgenii Kurkin) and O.L.; software, J.G.Q.-P. and E.K. (Evgenii Kurkin); validation, D.N. and J.G.Q.-P.; formal analysis, V.C., E.K. (Ekaterina Kurkina) and V.H.H.; resources, A.S.; writing—original draft preparation, A.S., J.G.Q.-P., O.L. and V.C.; writing—review and editing, J.G.Q.-P., E.K. (Evgenii Kurkin), O.L., D.N., V.C., E.K. (Ekaterina Kurkina), V.H.H. and A.S.; visualization, V.C., E.K. (Ekaterina Kurkina) and V.H.H.; supervision, D.N.; project administration, D.N.; funding acquisition, D.N., O.L. and A.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The raw data cannot be shared at this time as the data also forms part of an ongoing study.

Conflicts of Interest

The authors declare no conflicts of interest. The founders had no role in the design of the study; in the collection, analysis, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Acknowledgments

We would like to express our gratitude to Elena Tarasova and Damian Josue Guerra Guerra for the aerodynamic experiment assistance and Georgy Charkviani and Dmitry Kashirsky for help with propellers manufacturing.

References

1. G. Chen; D. Ma; Y. Jia; X. Xia; C. He Comprehensive Optimization of the Unmanned Tilt-Wing Cargo Aircraft with Distributed Propulsors., 2020, 8,pp. 137867-137883. DOI: https://doi.org/10.1109/ACCESS.2020.3012481.

2. L.W. Traub Considerations in optimal propeller design., 2021, 58,p. 8. DOI: https://doi.org/10.2514/1.C036258.

3. A. Betz, Pergamon Press: Oxford, UK, 1966,

4. A. Betz, Akademie der Wissenschaften: Gottingen, Germany, 1919,pp. 193-217.

5. C.N. Adkins; R.H. Liebeck Design of Optimum Propellers., 1994, 10,pp. 676-682. DOI: https://doi.org/10.2514/3.23779.

6. V.L. Aleksandrov, Gosudarstvennoe Izdatel’stvo Oboronnoj Promyshlennosti: Moscow, Russia, 1951,. (In Russian)

7. A. Kümmel; C. Breitsamter Multi-Disciplinary Framework for Propeller Blade Design., 2021, 1024,p. 012060. DOI: https://doi.org/10.1088/1757-899X/1024/1/012060.

8. W. Chen; K. Chiu; M.D. Fuge Airfoil Design Parameterization and Optimization Using Bézier Generative Adversarial Networks., 2020, 58,pp. 4723-4735. DOI: https://doi.org/10.2514/1.J059317.

9. O. Gur; A. Rosen Optimization of Propeller Bases Propulsion Systems., 2009, 46,pp. 95-106. DOI: https://doi.org/10.2514/1.36055.

10. J. Dorfling; K. Rokhsaz Constrained and Unconstrained Propeller Blade Optimization., 2015, 52,pp. 1179-1188. DOI: https://doi.org/10.2514/1.C032859.

11. P. Yu; J. Peng; J. Bai; X. Han; X. Song Aeroacoustic and Aerodynamic Optimization of Propeller Blades., 2020, 33,pp. 826-839. DOI: https://doi.org/10.1016/j.cja.2019.11.005.

12. S. Wang; S. Zhang; S. Ma An Energy Efficiency Optimization Method for Fixed Pitch Propeller Electric Aircraft Propulsion Systems., 2019, 7,pp. 159986-159993. DOI: https://doi.org/10.1109/ACCESS.2019.2950453.

13. E.G. Bekele; J.W. Nicklow Multi-objective automatic calibration of SWAT using NSGA-II., 2007, 341,pp. 165-176. DOI: https://doi.org/10.1016/j.jhydrol.2007.05.014.

14. R. Ma; B.W. Zhong; P.Q. Liu Optimization design study of low-Reynolds-number high-lift airfoils for the high-efficiency propeller of low-dynamic vehicles in stratosphere., 2010, 53,pp. 2792-2807. DOI: https://doi.org/10.1007/s11431-010-4087-0.

15. S. Slavik; J. Klesa; J. Brabec Propeller Selection by Means of Pareto-Optimal Sets Applied to Flight Performance., 2020, 7, 21. DOI: https://doi.org/10.3390/aerospace7030021.

16. B.R. Fang, Chinese Aviation Industry Press: Beijing, China, 1997,pp. 97-130. DOI: https://doi.org/10.3390/aerospace9020079.

17. A. Colozza, Federal Data Systems: Cleveland, OH, USA, 1998,. NASA/CR 98-208520

18. J. Morgado; M. Abdollahzadeh; M. Silvestre; J. Páscoa High Altitude Propeller Design and Analysis., 2015, 45,pp. 398-407. DOI: https://doi.org/10.1016/j.ast.2015.06.011.

19. T.A. Burdett; K.W. Van Treuren A Theoretical and Experimental Comparison of Optimizing Angle of Twist Using BET and BEMT.,pp. 797-809. DOI: https://doi.org/10.1115/GT2012-68350.

20. H.A. Oliveira; J.G. de Matos; L.A.d.S. Ribeiro; O.R. Saavedra; J.R.P. Vaz Assessment of Correction Methods Applied to BEMT for Predicting Performance of Horizontal-Axis Wind Turbines., 2023, 15, 7021. DOI: https://doi.org/10.3390/su15087021.

21. S. Goldstein On The Vortex Theory of Screw Propellers., 1929, 123,pp. 440-465. DOI: https://doi.org/10.1098/rspa.1929.0078.

22. W. Zhong; T.G. Wang; W.J. Zhu; W.Z. Shen Evaluation of Tip Loss Corrections to AD/NS Simulations of Wind Turbine Aerodynamic Performance., 2019, 9, 4919. DOI: https://doi.org/10.3390/app9224919.

23. E. Branlard; K. Dixon; M. Gaunaa Vortex Methods to Answer The Need For Improved Understanding And Modelling of Tip-Loss Factors., 2013, 7,pp. 311-320. DOI: https://doi.org/10.1049/iet-rpg.2012.0283.

24. Y. Elkhchine; M. Sriti Tip Loss Factor Effects on Aerodynamic Performances of Horizontal Axis Wind Turbine., 2017, 118,pp. 136-140. DOI: https://doi.org/10.1016/j.egypro.2017.07.028.

25. V.I. Shaidakov, MAI: Moscow, Russia, 1996,. (In Russian)

26. Q.R. Wald The Aerodynamics of Propellers., 2006, 42,pp. 85-128. DOI: https://doi.org/10.1016/j.paerosci.2006.04.001.

27. XFOIL Subsonic Airfoil Development System.. Available online: https://web.mit.edu/drela/Public/web/xfoil/ <date-in-citation content-type="access-date" iso-8601-date="2024-01-10">(accessed on 10 January 2024)</date-in-citation>.

28. J. Munguia; W. Van Treuren, Baylor University: Waco, TX, USA, 2019,

29. S.S. Pierret; R.A. Van den Braembussche Turbomachinery Blade Design Using a Navier-Stokes Solver and Artificial Neural Net-work., 1999, 121,pp. 326-332. DOI: https://doi.org/10.1115/1.2841318.

30. M.M. Segui; Y. Castelar; R.M. Botez Wing Airfoils Generation Based on a New Parametric Curve for Aerodynamic Optimization Application., Canadian Aeronautics and Space Institute: Montreal, QC, Canada, 2019,

31. A. Shikhar Jaiswal Shape Parameterization of Airfoil Shapes Using Bezier Curves., Springer: Berlin/Heidelberg, Germany, 2017, DOI: https://doi.org/10.1007/978-981-10-1771-1_13.. Lecture Notes in Mechanical Engineering

32. R.A. Alexeev; V.A. Tishchenko; V.G. Gribin; I.Y. Gavrilov Turbine Blade Profile Design Method Based on Bezier Curves., 2017, 891,p. 012254. DOI: https://doi.org/10.1088/1742-6596/891/1/012254.

33. X. Hao; W. Zhang; X. Liu; J. Liu Aerodynamic and Aeroacoustic Optimization of Wind Turbine Blade by a Genetic Algorithm., DOI: https://doi.org/10.2514/6.2008-1331.

34. A.I. Borovkov; I.B. Voinov; D.F. Ibraev Determination of the Optimal Aerodynamic Shape for a Propeller Blade Based on Parametric Optimization., 2021, 64,pp. 173-180. DOI: https://doi.org/10.3103/S106879982102001X.

35. H.C. Ayaz Ümütlü; Z. Kiral Airfoil Shape Optimization using Bézier Curve and Genetic Algorithm., 2022, 26,pp. 32-40. DOI: https://doi.org/10.3846/aviation.2022.16471.

36. P. Xin; L. Dawei; S. Jixiang; L. Yonghong Airfoil Aerodynamic Optimization Based on an Improved Genetic Algorithm.,pp. 133-137. DOI: https://doi.org/10.1109/ISDEA.2014.37.

37. J. Jiao; B.-F. Song; Y.-G. Zhang; Y.-B. Li Optimal Design and Experiment of Propellers for High Altitude Airship., 2017, 232,pp. 1887-1902. DOI: https://doi.org/10.1177/0954410017704217.

38. L.S. Bocii; L.P. Di Noia; R. Rizzo Optimization of the Energy Storage of Series-Hybrid Propelled Aircraft by means of Integer Differential Evolution., 2019, 6, 59. DOI: https://doi.org/10.3390/aerospace6050059.

39. A. Muratoglu; R. Tekin; O.F. Ertugrul Hydrodynamic Optimization of High-Performance Blade Sections for Stall Regulated Hydrokinetic Turbines Using Differential Evolution Algorithm., 2021, 220,p. 108389. DOI: https://doi.org/10.1016/j.oceaneng.2020.108389.

40. D.A. Spera, GRC: Cleveland, OH, USA, 2008,. Technical Report NASA/CR-2008-215434

41. J.G. Quijada Pioquinto; E.I. Kurkin; D.V. Nazarov; O.E. Lukyanov; V.O. Chertykovtseva, RU 2024610972, FIPS: Moscow, Russia, 2024,. Available online: https://new.fips.ru/registers-doc-view/fips_servlet?DB=EVM&DocNumber=2024610972&TypeFile=html <date-in-citation content-type="access-date" iso-8601-date="2024-02-20">(accessed on 20 February 2024)</date-in-citation>.

42. M.M. Ali; W.X. Zhu A penalty function-based differential evolution algorithm for constrained global optimization., 2013, 54,pp. 707-739. DOI: https://doi.org/10.1007/s10589-012-9498-3.

43. T. Van To; H. Ngoc Phien Development of Bezier-based curves., 1992, 20,pp. 109-115. DOI: https://doi.org/10.1016/0166-3615(92)90132-7.

44. V.V. Zherejov; A.N. Kusumov, KGTU: Kazan, Russia, 1997,. (In Russian)

45. A.S. Kravec, Uchebnoe Posobie, Gos. Izd. Oboronnoj Promyshlennosti: Moscow, Russia, 1941,. (In Russian)

46. R.W. Derksen; T. Rogalsky Bezier-PARSEC: An Optimized Aerofoil Parameterization for Design., 2010, 41,pp. 923-930. DOI: https://doi.org/10.1016/j.advengsoft.2010.05.002.

47. O.U. Espinosa Barcenas; J.G. Quijada Pioquinto; E. Kurkina; O. Lukyanov Surrogate aerodynamic wing modeling based on a multilayer perceptron., 2023, 10, 149. DOI: https://doi.org/10.3390/aerospace10020149.

48. J. Lin; H.D. Yao; C. Wang; Y. Su; C. Yang Hydrodynamic performance of a rim-driven thruster improved with gap geometry adjustment., 2023, 17,p. 2183902. DOI: https://doi.org/10.1080/19942060.2023.2183902.

49. T. Eltaeib; A. Mahmood Differential Evolution: A Survey and Analysis., 2018, 8, 1954. DOI: https://doi.org/10.3390/app8101945.

50. R. Tanabe; A.S. f*ckunaga Improving the search performance of SHADE using linear population size reduction.,pp. 1200-1207. DOI: https://doi.org/10.1109/CEC.2014.6900380.

51. O.U. Espinosa Barcenas; J.G. Quijada Pioquinto; E. Kurkina; O. Lukyanov Multidisciplinary analysis and optimization method for conceptually designing of electric flying-wing unmanned aerial vehicles., 2022, 6, 307. DOI: https://doi.org/10.3390/drones6100307.

52. I. Wong; W. Liu; C.H. Ho; X. Ding Continuous adaptive population reduction (CAPR) for differential evolution optimization., 2017, 22,pp. 289-305. DOI: https://doi.org/10.1177/2472630317690318. PMID: https://www.ncbi.nlm.nih.gov/pubmed/28378610.

53. R.L. Iman, John Wiley Sons: Hoboken, NJ, USA, 2008,pp. 408-411.

54. K. Zielinksi; P. Weitkemper; R. Laur; K.D. Kammeyer Examination of Stopping Criteria for Differential Evolution Based on a Power Allocation Problem.. Available online: https://api.semanticscholar.org/CorpusID:13358798 <date-in-citation content-type="access-date" iso-8601-date="2023-11-02">(accessed on 2 November 2023)</date-in-citation>.

55. H. Sharma Lightweight Pipelining in Python. Using Joblib for Storing the Machine Learning Pipeline to a File.. Available online: https://towardsdatascience.com/lightweight-pipelining-in-python-1c7a874794f4 <date-in-citation content-type="access-date" iso-8601-date="2023-11-02">(accessed on 2 November 2023)</date-in-citation>.

56. A.P. Piotrowski Review of differential evolution population size., 2017, 32,pp. 1-24. DOI: https://doi.org/10.1016/j.swevo.2016.05.003.

57. O.E. Lukyanov; O.U. Espinosa Barcenas; D.V. Zolotov Experimental model of an electric power plant for small UAV’s automatic control systems.,

Figures and Tables

Figure 1: The Bezier curves determine the variation of the design variables. [Please download the PDF to view the image]

Figure 2: Airfoil construction using cubic Bezier curves. [Please download the PDF to view the image]

Figure 3: Drawing of a blade from Bezier curves and the geometric twist curve. [Please download the PDF to view the image]

Figure 4: Aerodynamic coefficients of the CLARK Y profile by different methods. [Please download the PDF to view the image]

Figure 5: Calculation and interpolation of aerodynamic characteristics of the profile by propeller span. [Please download the PDF to view the image]

Figure 6: NACA 5868-9 propeller calculation in ANSYS CFX, (a) the problem statement, (b) the domain dimensions. [Please download the PDF to view the image]

Figure 7: Example of NACA 5868-9 propeller calculation for 1500 rpm and f[sub.0.75] = 25° case, (a) velocity in stationary domain, (b) pressure at propeller wall. [Please download the PDF to view the image]

Figure 8: Comparison by thrust coefficient (a) and power coefficient (b). [Please download the PDF to view the image]

Figure 9: Convergence of the metrics in the tests corresponding to case 1. [Please download the PDF to view the image]

Figure 10: Top view of the propeller blades obtained for case 1. [Please download the PDF to view the image]

Figure 11: Examples of the cross-sections of the blades obtained for case 1. The airfoils are normalized with chord of length 1. [Please download the PDF to view the image]

Figure 12: Convergence of the metrics in the tests corresponding to case 2. [Please download the PDF to view the image]

Figure 13: Top view of the propeller blades obtained for case 2. [Please download the PDF to view the image]

Figure 14: Examples of the cross-sections (airfoils) of the blades obtained for case 2. The airfoils are normalized with chord of length 1. [Please download the PDF to view the image]

Figure 15: (a) Injection mold for the manufacture of the propellers; (b) manufactured propellers. [Please download the PDF to view the image]

Figure 16: Experimental setup. [Please download the PDF to view the image]

Figure 17: The designed (a) and commercial (b) propellers installed in the experimental setup. [Please download the PDF to view the image]

Figure 18: Experimental data of thrust from required power dependence for propellers: (a) designed, (b) commercial DJI propeller. [Please download the PDF to view the image]

Figure 19: Mean measurement values and scatter of propeller thrust from power requirement: (a) comparison for designed propeller between experimental values and numerical test, calculated by ISM, (b) comparison between experimental values for designed and commercial DJI propeller. [Please download the PDF to view the image]

Table 1: Intervals of the design variables.

ParameterIntervalVariableInterval
Case 1Case 2Case 1Case 2

(c/d)[sub.r]

[0.03, 0.10]

[0.05, 0.07]

x[sub.tt]

[0.30, 0.40]

[0.30, 0.40]

(c/d)[sub.m]

[0.03, 0.10]

[0.08, 0.13]

r¯ [sub.xtm]

[0.30, 0.80]

[0.20, 0.50]

(c/d)[sub.t]

[0.01, 0.02]

[0.01, 0.03]

y[sub.cr]

[0.01, 0.05]

[0.05, 0.08]

r¯ [sub.cm]

[0.35, 0.60]

[0.20, 0.50]

y[sub.cm]

[0.01, 0.05]

[0.05, 0.08]

a[sub.r] [°]

[-7, 7]

[0, 5]

y[sub.ct]

[0.005, 0.05]

[0.05, 0.08]

a[sub.m] [°]

[-5, 7]

[0, 5]

r¯ [sub.ycm]

[0.30, 0.80]

[0.20, 0.50]

a[sub.t] [°]

[-5, 7]

[0, 5]

x[sub.cr]

[0.30, 0.40]

[0.30, 0.40]

r¯ [sub.am]

[0.25, 0.75]

[0.20, 0.50]

x[sub.cm]

[0.30, 0.45]

[0.30, 0.40]

y[sub.tr]

[0.12, 0.20]

[0.10, 0.20]

x[sub.ct]

[0.30, 0.45]

[0.30, 0.40]

y[sub.tm]

[0.12, 0.16]

[0.08, 0.10]

r¯ [sub.xcm]

[0.30, 0.80]

[0.20, 0.50]

y[sub.tt]

[0.10, 0.12]

[0.08, 0.10]

n[sub.m] [rev/min]

[5 × 10[sup.3], 10 × 10[sup.3]]

[5 × 10[sup.3], 10 × 10[sup.3]]

r¯ [sub.ytm]

[0.30, 0.80]

[0.20, 0.50]

B

[2, 4]

[2, 3]

x[sub.tr]

[0.30, 0.40]

[0.30, 0.40]

d [m]

[0.3, 0.3]

[0.254, 0.254]

x[sub.tm]

[0.30, 0.45]

[0.30, 0.40]

Table 2: Results obtained from case 1.

ParameterNP = 50NP = 100
Test 1Test 2Test 3MeanCV, [%]Test 4Test 5Test 6MeanCV, [%]

(c/d)[sub.r]

0.043

0.082

0.094

0.073

29.45

0.035

0.038

0.043

0.039

8.04

(c/d)[sub.m]

0.100

0.094

0.088

0.094

5.21

0.085

0.087

0.100

0.091

7.40

(c/d)[sub.t]

0.012

0.015

0.018

0.015

15.33

0.020

0.020

0.019

0.020

3.42

r¯ [sub.cm]

0.509

0.400

0.600

0.503

16.29

0.551

0.587

0.522

0.553

4.82

a[sub.r] [°]

0.243

-0.784

6.288

1.916

162.84

4.214

-7.000

6.245

1.153

505.15

a[sub.m] [°]

6.144

5.278

4.307

5.243

14.31

6.956

6.081

5.609

6.216

8.98

a[sub.t] [°]

4.823

6.973

2.324

4.707

40.36

4.632

4.770

5.217

4.873

5.12

r¯ [sub.am]

0.253

0.251

0.250

0.251

0.60

0.253

0.286

0.605

0.381

41.62

y[sub.tr]

0.070

0.060

0.066

0.066

6.41

0.060

0.061

0.061

0.061

0.75

y[sub.tm]

0.060

0.061

0.061

0.061

0.66

0.060

0.060

0.060

0.060

0.00

y[sub.tt]

0.059

0.060

0.054

0.057

3.99

0.053

0.060

0.051

0.055

6.71

r¯ [sub.ytm]

0.745

0.300

0.300

0.448

46.74

0.598

0.530

0.364

0.497

19.73

x[sub.tr]

0.327

0.300

0.309

0.312

3.52

0.333

0.320

0.308

0.320

3.21

x[sub.tm]

0.329

0.335

0.300

0.321

4.73

0.332

0.336

0.315

0.328

2.73

x[sub.tt]

0.330

0.315

0.306

0.317

3.15

0.339

0.384

0.345

0.356

5.65

r¯ [sub.xtm]

0.787

0.722

0.300

0.603

35.75

0.300

0.300

0.357

0.319

8.35

y[sub.cr]

0.050

0.010

0.029

0.030

54.70

0.026

0.047

0.011

0.028

52.36

y[sub.cm]

0.010

0.011

0.014

0.012

15.13

0.010

0.010

0.010

0.010

0.00

y[sub.ct]

0.005

0.016

0.028

0.016

57.93

0.006

0.005

0.006

0.006

8.31

r¯ [sub.ycm]

0.358

0.459

0.652

0.490

24.90

0.483

0.358

0.712

0.518

28.33

x[sub.cr]

0.338

0.399

0.331

0.356

8.65

0.398

0.395

0.364

0.386

3.95

x[sub.cm]

0.443

0.443

0.306

0.397

16.28

0.352

0.386

0.387

0.37

4.33

x[sub.ct]

0.361

0.441

0.310

0.370

14.59

0.399

0.367

0.422

0.396

5.77

r¯ [sub.xcm]

0.692

0.407

0.778

0.626

25.32

0.643

0.800

0.505

0.649

18.58

n[sub.m] [rev/min]

6156

6238

6765

6386

4.22

6281

6312

5993

6195

2.32

B

2

2

2

2

2

2

2

2

d [m]

0.300

0.300

0.300

0.300

0.300

0.300

0.300

0.300

W [W]

226.8

225.7

228.1

226.9

0.44

226.0

226.7

225.7

226.1

0.18

T [N]

7.513

7.504

7.506

7.508

0.05

7.501

7.512

7.512

7.509

0.07

?[sub.d]

0.828

0.831

0.823

0.827

0.41

0.830

0.828

0.832

0.830

0.18

Computation time [s]

3897

4728

4417

4347

7.89

9987

9526

10,703

10,072

4.81

Table 3: Results obtained from the case 2.

ParameterTest 1Test 2Test 3MeanCV, [%]

(c/d)[sub.r]

0.063

0.070

0.055

0.063

9.78

(c/d)[sub.m]

0.130

0.120

0.125

0.125

3.27

(c/d)[sub.t]

0.030

0.028

0.030

0.029

3.21

r¯ [sub.cm]

0.456

0.500

0.478

0.478

3.76

a[sub.r] [°]

0.000

0.718

3.579

1.432

107.93

a[sub.m] [°]

4.625

2.254

3.772

3.550

27.62

a[sub.t] [°]

1.619

1.193

3.734

2.182

50.92

r¯ [sub.am]

0.435

0.327

0.500

0.421

16.96

y[sub.tr]

0.063

0.073

0.093

0.076

16.34

y[sub.tm]

0.042

0.040

0.040

0.041

2.31

y[sub.tt]

0.046

0.043

0.042

0.044

3.89

r¯ [sub.ytm]

0.251

0.200

0.208

0.220

10.19

x[sub.tr]

0.300

0.382

0.356

0.346

9.89

x[sub.tm]

0.331

0.301

0.329

0.320

4.28

x[sub.tt]

0.361

0.300

0.390

0.350

10.71

r¯ [sub.xtm]

0.497

0.200

0.248

0.315

41.33

y[sub.cr]

0.064

0.062

0.069

0.065

4.53

y[sub.cm]

0.051

0.051

0.052

0.051

0.92

y[sub.ct]

0.050

0.053

0.051

0.051

2.43

r¯ [sub.ycm]

0.272

0.226

0.222

0.240

9.45

x[sub.cr]

0.386

0.398

0.334

0.372

7.45

x[sub.cm]

0.300

0.302

0.300

0.301

0.31

x[sub.ct]

0.323

0.301

0.301

0.308

3.36

r¯ [sub.xcm]

0.263

0.249

0.209

0.240

9.52

n[sub.m] [rev/min]

6720

7662

6705

7209

6.37

B

2

2

2

2

0.00

d [m]

0.254

0.254

0.254

0.254

0.00

W [W]

72.56

72.17

72.24

72.32

0.23

T [N]

6.505

6.516

6.505

6.509

0.08

?[sub.s]

0.649

0.654

0.652

0.652

0.32

Computation time [s]

3333

3743

3227

3434

6.48

Author Affiliation(s):

Institute of Aerospace Engineering, Samara National Research University, 34 Moskovskoe Shosse, 443086 Samara, Russia; [emailprotected] (E.K.); [emailprotected] (J.G.Q.-P.); [emailprotected] (O.L.); [emailprotected] (D.N.); [emailprotected] (V.C.); [emailprotected] (E.K.); [emailprotected] (V.H.H.)

Author Note(s):

[*] Correspondence: [emailprotected]

DOI: 10.3390/computation12030052

COPYRIGHT 2024 MDPI AG
No portion of this article can be reproduced without the express written permission from the copyright holder.

Copyright 2024 Gale, Cengage Learning. All rights reserved.


Algorithm for Propeller Optimization Based on Differential Evolution. (2024)
Top Articles
Latest Posts
Article information

Author: Kieth Sipes

Last Updated:

Views: 5907

Rating: 4.7 / 5 (67 voted)

Reviews: 82% of readers found this page helpful

Author information

Name: Kieth Sipes

Birthday: 2001-04-14

Address: Suite 492 62479 Champlin Loop, South Catrice, MS 57271

Phone: +9663362133320

Job: District Sales Analyst

Hobby: Digital arts, Dance, Ghost hunting, Worldbuilding, Kayaking, Table tennis, 3D printing

Introduction: My name is Kieth Sipes, I am a zany, rich, courageous, powerful, faithful, jolly, excited person who loves writing and wants to share my knowledge and understanding with you.