Algorithm for Propeller Optimization Based on Differential Evolution

1. Introduction

The air propeller is one of the most important elements of an aircraft that provides lift and propulsion in the air. The relevance of using the propeller as a propulsion system has increased with the current growth in the production of unmanned aerial vehicles for various purposes around the world. An optimized air propeller can significantly reduce emissions, improve acoustic response, and enhance performance [1]. The design of air propellers to minimize propulsion system energy costs began with Zhukovsky, Betz, and Goldstein [2,3] and was then refined by a few scientists over the course of a century [4]. In recent years, computational fluid dynamics using CAE systems has become an important method for propeller design and computation [5,6,7]. Propeller optimization is a complex problem, often multi-objective optimization, requiring consideration of many factors and constraints, including not only aerodynamic issues but also strength, acoustics, etc. [8,9,10,11]. Therefore, the use of the finite volume method becomes irrational because it requires large computational power. 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 [12]. The paper [13] presents a methodology for inverse propeller blade design. The isolated blade section method is used to directly calculate propeller aerodynamic characteristics during the optimization process. This involved the calculation of the aerodynamic characteristics of the profiles in the blade sections using the xflr code. The obtained characteristics are verified with CFD calculations. The authors of [14] solved a multi-objective optimization problem subject to a set of constraints using a direct search algorithm. The propeller characteristics were calculated from the profiles in the blade sections using the discrete vortex method described in the Xfoil code. The paper [15] 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 propeller characteristics were calculated from the profiles in the blade sections using vortex methods. The optimization results showed that the desirable trade-off between propeller efficiency and weight is related to power consumption, structural strength, and even manufacturing and installation conditions. The authors [16] developed a procedure to design and fabricate propellers for small unmanned aerial vehicles. The method of sections implemented in QProp is the basis for the calculation of the direct problem of aerodynamics, and the optimization process is managed by the Bearcontrol 9 written in MATLAB, which allows to create a Visual Basic code for automatic construction of 3D model in SolidWorks.

The propeller has a complex three-dimensional geometry, and its optimization problem involves a parametric approach to define the geometric model. Bezier curves and surfaces are one of the most universal and robust approaches for specifying complex shapes [17] and can be used to determine the optimal shape of propellers [18]. This paper is devoted to the development of an efficient algorithm for the parametric optimization of an unmanned aerial vehicle propeller. This paper describes step-by-step the selection of methods of parametrization of geometric model, design calculation of propeller aerodynamics, and the use of metaheuristic algorithms of parametric optimization.

2. Materials and Methods

2.1. Mathematical Model Of Propeller Optimization

The objective of optimization is to determine the geometry and appropriate operating conditions of a propeller, which requires the least operating power, depending on the operating speed of the aircraft and the diameter of the propeller, also subject to a certain thrust. Mathematically, the optimization problem can be expressed as a constrained optimization problem:

min W(x)
such that T_min − 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_min is the minimum desired thrust; x is the vector of design parameters belonging to the feasible set of solutions X.

In order 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 [16]. The penalty function is expressed as:

$$\mathrm{L}\left(\mathrm{x}\right)=\left\{\begin{array}{cc}\mathrm{W}\left(\mathrm{x}\right)& \mathrm{i}\mathrm{f}\mathsf{\psi }\left(\mathrm{x}\right)=0\\ \mathrm{R}\mathsf{\psi }\left(\mathrm{x}\right)+\mathrm{U}\mathrm{*}& \mathrm{i}\mathrm{f}\mathsf{\psi }\left(\mathrm{x}\right)>0\wedge W\left(x\right)\le U*\\ \mathrm{R}\mathsf{\psi }\left(\mathrm{x}\right)+\mathrm{W}\left(\mathrm{x}\right)& \mathrm{i}\mathrm{f}\mathsf{\psi }\left(\mathrm{x}\right)>0\wedge W\left(x\right)>U*\end{array}\right.$$

(1)

where

$$\mathsf{\psi }\left(\mathrm{x}\right)=\max \left\{0,{\mathrm{T}}_{\mathrm{m}\mathrm{i}\mathrm{n}}-\mathrm{T}\left(\mathrm{x}\right)\right\}$$

(2)

where R is a penalty parameter; U* is an upper bound on the constrained global minimum value, which is provided initially by the user. ψ(x)>0 only if x is infeasible. 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_m, rev/min); and others that change depending on the radius of the propeller, such as the chord length, the effective angle of the airfoil and the geometry of the airfoil. It was considered that this type of variables has a continuous and continuous variation along the blade of the propeller. This variation was achieved by using two quadratic Bezier curves [19], (see Figure 1). A quadratic Bezier curve is the path drawn by the following function:

$$\mathrm{B}\left(\mathrm{t}\right)=\left(1-{\mathrm{t}}^2\right){\mathrm{P}}_0+2\left(1-\mathrm{t}\right)\mathrm{t}{\mathrm{P}}_1+{\mathrm{t}}^2{\mathrm{P}}_2,$$

(3)

where P₀, P₁ and P₂ are control points, and t is a parameter that always varies from 0 to 1.

B(t) can also be decomposed as (B_r̅(t), B_x(t)), where r̅ indicates the relative radius of the propeller and x the design parameter.

$${\mathrm{B}}_{\overline{\mathrm{r}}}\left(\mathrm{t}\right)={(1-\mathrm{t})}^2{\overline{\mathrm{r}}}_0+2\left(1-\mathrm{t}\right)\mathrm{t}{\overline{\mathrm{r}}}_1+{\mathrm{t}}^2{\overline{\mathrm{r}}}_2,$$

(4)

$${\mathrm{B}}_{\mathrm{x}}\left(\mathrm{t}\right)={(1-\mathrm{t})}^2{\mathrm{x}}_0+2\left(1-\mathrm{t}\right)\mathrm{t}{\mathrm{x}}_1+{\mathrm{t}}^2{\mathrm{x}}_2,$$

(5)

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

• Root curve

$$\left\{\begin{array}{c}{\overline{\mathrm{r}}}_0^{\mathrm{r}}=0.1\\ {\overline{\mathrm{r}}}_1^{\mathrm{r}}=0.5{\overline{\mathrm{r}}}_{\mathrm{x}\mathrm{m}}+0.05\\ {\overline{\mathrm{r}}}_2^{\mathrm{r}}={\overline{\mathrm{r}}}_{\mathrm{x}\mathrm{m}}\end{array}\right.\hspace{1em}\left\{\begin{array}{c}{\mathrm{x}}_0^{\mathrm{r}}={\mathrm{x}}_{\mathrm{r}}\\ {\mathrm{x}}_1^{\mathrm{r}}={\mathrm{x}}_{\mathrm{m}}\\ {\mathrm{x}}_2^{\mathrm{r}}={\mathrm{x}}_{\mathrm{m}}\end{array}\right.$$

(6)

• Tip curve

$$\left\{\begin{array}{c}{\overline{\mathrm{r}}}_0^{\mathrm{t}}={\overline{\mathrm{r}}}_{\mathrm{x}\mathrm{m}}\\ {\overline{\mathrm{r}}}_1^{\mathrm{t}}=0.5{\overline{\mathrm{r}}}_{\mathrm{x}\mathrm{m}}+0.485\\ {\overline{\mathrm{r}}}_2^{\mathrm{t}}=0.97\end{array}\right.\hspace{1em}\left\{\begin{array}{c}{\mathrm{x}}_0^{\mathrm{t}}={\mathrm{x}}_{\mathrm{m}}\\ {\mathrm{x}}_1^{\mathrm{t}}={\mathrm{x}}_{\mathrm{m}}\\ {\mathrm{x}}_2^{\mathrm{t}}={\mathrm{x}}_{\mathrm{t}}\end{array}\right.$$

(7)

In Equations 4 to 7, x may be the chord relative to the propeller diameter (c/d), the effective angle (α), maximum thickness (y_t), position of the maximum thickness relative to the chord (x_t), the maximum camber (y_c), and the position of the maximum camber relative to the chord (x_c). Each of these variables is with respect to the profile of each section of the propeller blade. The variable r̅_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_i = [(c/d)_r, (c/d)_m, (c/d)_t, r̅_cm, α_r, α_m, α_t, r̅_αm, y_tr, y_tm, y_tt, r̅_ytm, x_tr, x_tm, x_tt, r̅_xtm, y_cr, y_cm, y_ct, r̅_ycm, x_cr, x_cm, x_ct, r̅_xcm, n_m, B, d]_i.

2.3. Model For The Resolution Of The Objective Function

The objective function of the optimization process is the calculation of the required propeller power (W, in W) to achieve a thrust (T, in N). A modified version of the Isolated Sections Method (ISM) [20] was used to obtain these values. The variation between the original and our modified version of the ISM is to replace the geometric twisting of the sections (φ, in deg) by the effective angle of attack of the sections (α, in deg) as the input value.

The input variables required by the method are the following: B, V_∞, d, n(n, rev/s), and physical characteristics of the air (density (ρ, in kg/m³), kinematic viscosity (ν, [m²/s]), sound speed (a, [m/s])). ISM requires sectioning one of the propeller blades into an infinite 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_U, Y_U), (X_L, Y_L), each coordinate is in m).

To obtain the coordinates of the profile in each section, an airfoil parameterization method based on the Bezier-PARSEC technique was proposed [21]. Each airfoil is constructed with 4 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:

$${\mathrm{B}}_{\mathrm{x}}\left(\mathrm{t}\right)={(1-\mathrm{t})}^3{\mathrm{x}}_0+3{(1-\mathrm{t})}^2\mathrm{t}{\mathrm{x}}_1+3(1-\mathrm{t}){\mathrm{t}}^2{\mathrm{x}}_2+{\mathrm{t}}^3{\mathrm{x}}_3,$$

(8)

$${\mathrm{B}}_{\mathrm{y}}\left(\mathrm{t}\right)={(1-\mathrm{t})}^3{\mathrm{y}}_0+3{(1-\mathrm{t})}^2\mathrm{t}{\mathrm{y}}_1+3(1-\mathrm{t}){\mathrm{t}}^2{\mathrm{y}}_2+{\mathrm{t}}^3{\mathrm{y}}_3,$$

(9)

The construction of the profile is carried out using the following equations:

$${\mathrm{X}}_{\mathrm{t}}={\mathrm{B}}_{\mathrm{x}\mathrm{t}}^{\mathrm{l}\mathrm{e}}\left(\mathrm{t}\right)+{\mathrm{B}}_{\mathrm{x}\mathrm{t}}^{\mathrm{t}\mathrm{e}}\left(\mathrm{t}\right)$$

(10)

$${\mathrm{Y}}_{\mathrm{t}}={\mathrm{B}}_{\mathrm{y}\mathrm{t}}^{\mathrm{l}\mathrm{e}}\left(\mathrm{t}\right)+{\mathrm{B}}_{\mathrm{y}\mathrm{t}}^{\mathrm{t}\mathrm{e}}\left(\mathrm{t}\right)$$

(11)

$${\mathrm{X}}_{\mathrm{c}}={\mathrm{B}}_{\mathrm{x}\mathrm{c}}^{\mathrm{l}\mathrm{e}}\left(\mathrm{t}\right)+{\mathrm{B}}_{\mathrm{x}\mathrm{c}}^{\mathrm{t}\mathrm{e}}\left(\mathrm{t}\right)$$

(12)

$${\mathrm{Y}}_{\mathrm{c}}={\mathrm{B}}_{\mathrm{y}\mathrm{c}}^{\mathrm{l}\mathrm{e}}\left(\mathrm{t}\right)+{\mathrm{B}}_{\mathrm{y}\mathrm{c}}^{\mathrm{t}\mathrm{e}}\left(\mathrm{t}\right)$$

(13)

$$\mathsf{\theta }={\tan }^{-1}\left(\raisebox{1ex}{$\mathrm{d}{\mathrm{Y}}_{\mathrm{c}}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{d}{\mathrm{X}}_{\mathrm{c}}$}\right.\right)$$

(14)

The control points for the leading edge thickness curve are defined by:

$$\left\{\begin{array}{c}{\mathrm{x}}_0^{\mathrm{l}\mathrm{e}}=0\\ {\mathrm{x}}_1^{\mathrm{l}\mathrm{e}}=0\\ \begin{array}{c}{\mathrm{x}}_2^{\mathrm{l}\mathrm{e}}=0.5{\mathrm{x}}_{\mathrm{t}}\\ {\mathrm{x}}_3^{\mathrm{l}\mathrm{e}}={\mathrm{x}}_{\mathrm{t}}\end{array}\end{array}\right.\hspace{1em}\left\{\begin{array}{c}{\mathrm{y}}_0^{\mathrm{l}\mathrm{e}}=0\\ {\mathrm{y}}_1^{\mathrm{l}\mathrm{e}}=0.34{\mathrm{y}}_{\mathrm{t}}\\ \begin{array}{c}{\mathrm{y}}_2^{\mathrm{l}\mathrm{e}}={0.5\mathrm{y}}_{\mathrm{t}}\\ {\mathrm{y}}_3^{\mathrm{l}\mathrm{e}}=0.5{\mathrm{y}}_{\mathrm{t}}\end{array}\end{array}\right.$$

(15)

The control points for the trailing edge thickness curve are defined by:

$$\left\{\begin{array}{c}{\mathrm{x}}_0^{\mathrm{t}\mathrm{e}}={\mathrm{x}}_{\mathrm{t}}\\ {\mathrm{x}}_1^{\mathrm{t}\mathrm{e}}=0.3+0.7{\mathrm{x}}_{\mathrm{t}}\\ \begin{array}{c}{\mathrm{x}}_2^{\mathrm{t}\mathrm{e}}=0.6+0.4{\mathrm{x}}_{\mathrm{t}}\\ {\mathrm{x}}_3^{\mathrm{t}\mathrm{e}}=1\end{array}\end{array}\right.\hspace{1em}\left\{\begin{array}{c}{\mathrm{y}}_0^{\mathrm{t}\mathrm{e}}=0.5{\mathrm{y}}_{\mathrm{t}}\\ {\mathrm{y}}_1^{\mathrm{t}\mathrm{e}}=0.5{\mathrm{y}}_{\mathrm{t}}\\ \begin{array}{c}{\mathrm{y}}_2^{\mathrm{t}\mathrm{e}}={0.29\mathrm{y}}_{\mathrm{t}}\\ {\mathrm{y}}_3^{\mathrm{t}\mathrm{e}}=0\end{array}\end{array}\right.$$

(16)

The control points for the leading edge camber curve are defined by:

$$\left\{\begin{array}{c}{\mathrm{x}}_0^{\mathrm{l}\mathrm{e}}=0\\ {\mathrm{x}}_1^{\mathrm{l}\mathrm{e}}={\mathrm{x}}_{\mathrm{c}}/3\\ \begin{array}{c}{\mathrm{x}}_2^{\mathrm{l}\mathrm{e}}=2{\mathrm{x}}_{\mathrm{c}}/3\\ {\mathrm{x}}_3^{\mathrm{l}\mathrm{e}}={\mathrm{x}}_{\mathrm{c}}\end{array}\end{array}\right.\hspace{1em}\left\{\begin{array}{c}{\mathrm{y}}_0^{\mathrm{l}\mathrm{e}}=0\\ {\mathrm{y}}_1^{\mathrm{l}\mathrm{e}}=0.71{\mathrm{y}}_{\mathrm{c}}\\ \begin{array}{c}{\mathrm{y}}_2^{\mathrm{l}\mathrm{e}}={\mathrm{y}}_{\mathrm{c}}\\ {\mathrm{y}}_3^{\mathrm{l}\mathrm{e}}={\mathrm{y}}_{\mathrm{c}}\end{array}\end{array}\right.$$

(17)

And the control points for the trailing edge camber curve are defined by:

$$\left\{\begin{array}{c}{\mathrm{x}}_0^{\mathrm{t}\mathrm{e}}={\mathrm{x}}_{\mathrm{c}}\\ {\mathrm{x}}_1^{\mathrm{t}\mathrm{e}}=\left(1+2{\mathrm{x}}_{\mathrm{c}}\right)/3\\ \begin{array}{c}{\mathrm{x}}_2^{\mathrm{t}\mathrm{e}}=\left(2+{\mathrm{x}}_{\mathrm{c}}\right)/3\\ {\mathrm{x}}_3^{\mathrm{t}\mathrm{e}}=1\end{array}\end{array}\right.\hspace{1em}\left\{\begin{array}{c}{\mathrm{y}}_0^{\mathrm{t}\mathrm{e}}={\mathrm{y}}_{\mathrm{c}}\\ {\mathrm{y}}_1^{\mathrm{t}\mathrm{e}}={\mathrm{y}}_{\mathrm{c}}\\ \begin{array}{c}{\mathrm{y}}_2^{\mathrm{t}\mathrm{e}}={0.43\mathrm{y}}_{\mathrm{c}}\\ {\mathrm{y}}_3^{\mathrm{t}\mathrm{e}}=0\end{array}\end{array}\right.$$

(18)

The points for determining the upper curve of the airfoil are determined by:

$${\mathrm{X}}_{\mathrm{U}}={\mathrm{X}}_{\mathrm{c}}-{\mathrm{Y}}_{\mathrm{t}}\sin \mathsf{\theta }$$

(19)

$${\mathrm{Y}}_{\mathrm{U}}={\mathrm{Y}}_{\mathrm{c}}+{\mathrm{Y}}_{\mathrm{t}}\cos \mathsf{\theta }$$

(20)

And the points for determining the lower curve of the airfoil are determined by:

$${\mathrm{X}}_{\mathrm{L}}={\mathrm{X}}_{\mathrm{c}}+{\mathrm{Y}}_{\mathrm{t}}\sin \mathsf{\theta }$$

(21)

$${\mathrm{Y}}_{\mathrm{L}}={\mathrm{Y}}_{\mathrm{c}}-{\mathrm{Y}}_{\mathrm{t}}\cos \mathsf{\theta }$$

(22)

All the distribution curves of the design variables are generated by Algorithm 1.

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

$\mathbf{Inputs}:{\mathrm{x}}_{\mathrm{i}},\mathrm{d},\overline{\mathrm{r}}$, NS; Outputs: (c/d)i, αi, xti, yti, xci, yci, XU, YU, XL, YL; $$\begin{gathered} \mathrm{Create}\mathrm{the}\mathrm{distribution}\mathrm{curve}\mathrm{for}(\mathrm{c}/\mathrm{d}{)}_{\mathrm{i}}\mathrm{as}\mathrm{a}\mathrm{function}\mathrm{of}\mathrm{the}{\overline{\mathrm{r}}}_{\mathrm{i}}\mathrm{of}\mathrm{the}\mathrm{blade}\mathrm{with}(4),(5),(6),(7); \\ \mathrm{Create}\mathrm{the}\mathrm{distribution}\mathrm{curve}\mathrm{for}{\mathsf{\alpha }}_{\mathrm{i}}\mathrm{as}\mathrm{a}\mathrm{function}\mathrm{of}\mathrm{the}{\overline{\mathrm{r}}}_{\mathrm{i}}\mathrm{of}\mathrm{the}\mathrm{blade}\mathrm{with}(4),(5),(6),(7); \\ \mathrm{Create}\mathrm{the}\mathrm{distribution}\mathrm{curve}\mathrm{for}{\mathrm{x}}_{\mathrm{ti}}\mathrm{as}\mathrm{a}\mathrm{function}\mathrm{of}\mathrm{the}{\overline{\mathrm{r}}}_{\mathrm{i}}\mathrm{of}\mathrm{the}\mathrm{blade}\mathrm{with}(4),(5),(6),(7); \\ \mathrm{Create}\mathrm{the}\mathrm{distribution}\mathrm{curve}\mathrm{for}{\mathrm{y}}_{\mathrm{ti}}\mathrm{as}\mathrm{a}\mathrm{function}\mathrm{of}\mathrm{the}{\overline{\mathrm{r}}}_{\mathrm{i}}\mathrm{of}\mathrm{the}\mathrm{blade}\mathrm{with}(4),(5),(6),(7); \\ \mathrm{Create}\mathrm{the}\mathrm{distribution}\mathrm{curve}\mathrm{for}{\mathrm{x}}_{\mathrm{ci}}\mathrm{as}\mathrm{a}\mathrm{function}\mathrm{of}\mathrm{the}{\overline{\mathrm{r}}}_{\mathrm{i}}\mathrm{of}\mathrm{the}\mathrm{blade}\mathrm{with}(4),(5),(6),(7); \\ \mathrm{Create}\mathrm{the}\mathrm{distribution}\mathrm{curve}\mathrm{for}{\mathrm{y}}_{\mathrm{ci}}\mathrm{as}\mathrm{a}\mathrm{function}\mathrm{of}\mathrm{the}{\overline{\mathrm{r}}}_{\mathrm{i}}\mathrm{of}\mathrm{the}\mathrm{blade}\mathrm{with}(4),(5),(6),(7); \\ //\mathrm{Get}\mathrm{the}\mathrm{airfoil}\mathrm{in}\mathrm{each}\sec \mathrm{tion}\mathrm{of}\mathrm{the}\mathrm{blade} \\ \mathbf{for}\mathrm{s}=1\mathbf{to}\mathrm{NS}\mathbf{do} \end{gathered}$$       $\mathrm{Get}{\mathrm{x}}_{\mathrm{t}}\left({\overline{\mathrm{r}}}_{\mathrm{s},\mathrm{i}}\right),{\mathrm{y}}_{\mathrm{t}}\left({\overline{\mathrm{r}}}_{\mathrm{s},\mathrm{i}}\right),{\mathrm{x}}_{\mathrm{c}}\left({\overline{\mathrm{r}}}_{\mathrm{s},\mathrm{i}}\right)\mathrm{and}{\mathrm{y}}_{\mathrm{c}}\left({\overline{\mathrm{r}}}_{\mathrm{s},\mathrm{i}}\right);$    Get Xts,i and Yts,i with (8), (9), (10), (11), (15), (16);    Get Xcs,i and Ycs,i with (8), (9), (12), (13), (17), (18);    Get θs,i with (14);    Create and save the points of the sth-airfoil XU, YU, XL, YL with (19), (20), (21), (22);

The iterative method to obtain the thrust provided by the propeller and the required power is shown in [20]. This method was modified to use the effective angle of attack (α) of the profile in each section as a constant and the geometric twist (φ) of the section as a variable that is updated in each iteration.

To obtain the aerodynamic coefficients of each section of the blade, it was obtained by using the Xfoil program, which shows a good performance for the evaluation of profiles (see Figure 3).

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.

$$\mathsf{\omega }=2\mathsf{\pi }\mathrm{n}$$

(23)

$$\overline{\mathrm{v}}=\frac{{\mathrm{V}}_{\mathrm{\infty }}}{\mathsf{\omega }\mathrm{R}}$$

(24)

The 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 α of each section as input, it is necessary to calculate in the first instance 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 to reduce the number of calculations, interpolation of the profile data was performed over a smaller number of sections, in each of which the characteristics of the angle of attack were calculated considering the local values of the 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.

Having the local characteristics of the flow, in addition to the chord, geometry and angle of attack of the airfoil, coefficients of lift (c_l) and drag (c_d) coefficients in each section are calculated. The next step is to get ${\overline{\mathrm{U}}}_1$ , ${\overline{\mathrm{V}}}_1$ and ${\overline{\mathsf{\Gamma }}}_{\mathrm{r}}$ by an iterative process, which are necessary to get the coefficients ct and mk. The equations involved in this process are shown below:

$${\overline{\mathrm{v}}}_1=-\frac{\overline{\mathrm{v}}}{2}+\sqrt{\frac{{\overline{\mathrm{v}}}^2}{4}+{\overline{\mathrm{u}}}_1\left(\overline{\mathrm{r}}-{\overline{\mathrm{u}}}_1\right)+2\underset{\overline{\mathrm{r}}}{\overset{1}{\int }}\frac{{\overline{\mathrm{u}}}_1^2}{\overline{\mathrm{r}}}\mathrm{d}\overline{\mathrm{r}}}$$

(25)

$${\overline{\mathrm{U}}}_1=\overline{\mathrm{r}}-{\overline{\mathrm{u}}}_1$$

(26)

$${\overline{\mathrm{V}}}_1=\overline{\mathrm{v}}+{\overline{\mathrm{v}}}_1$$

(27)

$${\overline{\mathrm{W}}}_1=\sqrt{{\overline{\mathrm{V}}}_1^2+{\overline{\mathrm{U}}}_1^2}$$

(28)

$${\mathsf{\beta }}_1={\tan }^{-1}\left(\frac{{\overline{\mathrm{V}}}_1}{{\overline{\mathrm{U}}}_1}\right)$$

(29)

$$\mathsf{\sigma }=\frac{\mathrm{B}\mathrm{c}}{\mathrm{R}\mathsf{\pi }}$$

(30)

$${\overline{\mathsf{\Gamma }}}_{\mathrm{r}}=\frac{1}{8}\mathsf{\sigma }{\mathrm{c}}_{\mathrm{l}}{\overline{\mathrm{W}}}_1$$

(31)

$${\mathrm{f}}_{\mathrm{r}}=\frac{2}{\mathsf{\pi }}{\cos }^{-1}\left({\mathrm{e}}^{-\frac{0.5\mathrm{B}(1-\overline{\mathrm{r}})}{\overline{\mathrm{r}}\sin \left({\mathsf{\beta }}_1\right)}}\right)$$

(32)

previously before initializing the iterative process, it is important to initialize the n values of ${\overline{\mathrm{u}}}_1$ and ${\int }_{\overline{r}}^1\frac{{\overline{u}}_1^2}{\overline{r}}d\overline{r}$ equal to zero. At the end of each iteration ${\overline{\mathrm{u}}}_1$ has to be updated making use by:

$${\overline{\mathrm{u}}}_1=\frac{{\overline{\mathsf{\Gamma }}}_{\mathrm{r}}}{{\mathrm{f}}_{\mathrm{r}}\overline{\mathrm{r}}}$$

(33)

while the values of ${\int }_{\overline{r}}^1\frac{{\overline{u}}_1^2}{\overline{r}}d\overline{r}$ 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_t and m_k, it is first necessary to calculate their differentials in each section, and then integrate with respect to r̅, making use of trapz().

$${\mathrm{d}\mathrm{c}}_{\mathrm{t}}=8{\overline{\mathsf{\Gamma }}}_{\mathrm{r}}\left({\overline{\mathrm{U}}}_1-{\mathrm{K}}^{-1}{\overline{\mathrm{V}}}_1\right)$$

(34)

$${\mathrm{d}\mathrm{m}}_{\mathrm{k}}=8{\overline{\mathsf{\Gamma }}}_{\mathrm{r}}\left({\overline{\mathrm{V}}}_1+{\mathrm{K}}^{-1}{\overline{\mathrm{U}}}_1\right)\overline{\mathrm{r}}$$

(35)

By obtaining the coefficients c_t and m_k, the thrust provided by the propeller (T) and the required power of the propeller (W) can be calculated.

$$\mathrm{T}=0.5{\mathrm{c}}_{\mathrm{t}}\mathsf{\rho }{\left(\mathsf{\omega }\mathrm{R}\right)}^2\mathsf{\pi }{\mathrm{R}}^2$$

(36)

$$\mathrm{W}=0.5{\mathrm{m}}_{\mathrm{k}}\mathsf{\rho }{\left(\mathsf{\omega }\mathrm{R}\right)}^3\mathsf{\pi }{\mathrm{R}}^2$$

(37)

Other propeller performance metrics that can be obtained with this method are the dynamic efficiency ${\mathsf{\eta }}_{\mathrm{d}}$ and the static efficiency ${\mathsf{\eta }}_{\mathrm{s}}$ of the propeller.

$${\mathsf{\alpha }}_{\mathrm{p}}=\frac{\mathrm{T}}{\mathsf{\rho }{\mathrm{n}}^2{\mathrm{d}}^4}$$

(38)

$${\mathsf{\beta }}_{\mathrm{p}}=\frac{\mathrm{W}}{\mathsf{\rho }{\mathrm{n}}^3{\mathrm{d}}^5}$$

(39)

$${\mathsf{\lambda }}_{\mathrm{p}}=\frac{{\mathrm{V}}_{\mathrm{\infty }}}{\mathrm{n}\mathrm{d}}$$

(40)

$${\mathsf{\eta }}_{\mathrm{d}}=\frac{{\mathsf{\alpha }}_{\mathrm{p}}{\mathsf{\lambda }}_{\mathrm{p}}}{{\mathsf{\beta }}_{\mathrm{p}}}$$

(41)

$${\mathsf{\eta }}_{\mathrm{s}}=\frac{{\mathrm{c}}_{\mathrm{t}}^{3/2}}{2{\mathrm{m}}_{\mathrm{k}}}$$

(42)

Algorithm 2: Modified ISM
Input n, B, V∞, d, ρ, ν, a, c, $\overline{r}$, [XU, YU, XL, YL], NS;
R = d/2;
Get ω with (33);
Get $\overline{v}$ with (34);
for s = 1 to NS do
	rs = ${\overline{r}}_s$R;
	Get Res and Ms in each section;
	cls, cds = runXFOIL(αs, [XU, YU, XL, YL]s, Res, Ms);
	0 →${\overline{u}}_{1_s}$, 0→${\left({\int }_{\overline{r}}^1\frac{{\overline{u}}_1^2}{\overline{r}}d\overline{r}\right)}_s$
for l = 1 to 10 do
	Iu = ⌀;;
	for s = 1 to NS do
		Get ${\overline{v}}_{1_s}$ with (25);
		Get ${\overline{U}}_{1_s}$ with (26);
		Get ${\overline{V}}_{1_s}$ with (27);
		Get ${\overline{W}}_{1_s}$ with (28);
		Get β1s with (29);
		Φs = αs + β1s;
		Ks = cls/cds;
		Get σs with (30);
		Get ${\overline{\Gamma }}_{r_s}$ with (31);
		Get frs with (32);
		Update ${\overline{u}}_{1_s}$ with (33);
		${\left(\frac{{\overline{u}}_{1_s}^2}{{\overline{r}}_s}\right)}_s$→$I_{u_s}$;
	fors = 1 to NS do
		Get ${\left({\int }_{\overline{r}}^1\frac{{\overline{u}}_1^2}{\overline{r}}d\overline{r}\right)}_s$ with trapz(Iu[s:end], $\overline{r}$[s:end]);
fors = 1 to NS do
	Get ${dc}_{t_s}$ with (34);
	Get ${dm}_{k_s}$ with (35);
ct = trapz(dct, $\overline{r}$);
mk = trapz(dmk, $\overline{r}$);
Get αp, βp and λp with (38), (39) and (40) respectively;
Get outputs T, W, ηd, and ηs with (36), (37), (41) and (42) respectively;

Verification of the proposed method of propeller calculation has been 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 [22]), as well as by comparison with the results of numerical mathematical modeling by solving the Navier-Stokes equations in CFX software. The results of comparison in terms of thrust coefficient and power are shown in Figure 5.

2.4. 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 the evolutionary operators, methods of population size reduction (PSR), sampling techniques for 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.

DE 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) [23]. As with other evolutionary algorithms, the search performance of DE algorithms depends on control parameter settings. A standard DE has three main control parameters, which are the 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 DE 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 the online control parameters during the search process. One of the variants of DE that apply this type of schemes is the success-history based adaptation for DE (SHADE) [24].

SHADE uses a historical memory MCR, MF 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.

Algorithm 3: Memory update algorithm in SHADE
IfSCR ≠ ∅ and SF ≠ ∅ then
	If MCR,k,g = -1 ormax(SCR) = 0 then
		MCR,k,g+1 = -1;
	else
		MCR,k,g+1 = meanWL(SCR);
	MF,k,g+1 = meanWL(SF);
	k++;
	If k > Hthen, k = 1;
else
	MCR,k,g+1 = MCR,k,g;
	MF,k,g+1 = MF,k,g;

In Algorithm 3, index k (1 ≤ k ≤ H) determines the position in the memory to update. In generation g, the k-th 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 update algorithm 1, note that when all individuals in generation g fail to generate a trial vector which is better than the parent, i.e., SCR = SF = ∅, the memory is not updated. The weighted Lehmer mean meanWL(S) is computed using the formula below:

$${\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}_{\mathrm{W}\mathrm{L}}\left(\mathrm{S}\right)=\frac{{\sum }_{\mathrm{m}=1}^{\left|\mathrm{S}\right|}{\mathrm{w}}_{\mathrm{m}}{\mathrm{S}}_{\mathrm{m}}^2}{{\sum }_{\mathrm{m}=1}^{\left|\mathrm{S}\right|}{\mathrm{w}}_{\mathrm{m}}{\mathrm{S}}_{\mathrm{m}}}$$

(43)

$${\mathrm{w}}_{\mathrm{m}}=\frac{\mathsf{\Delta }{\mathrm{f}}_{\mathrm{m}}}{{\sum }_{\mathrm{l}=1}^{\left|\mathrm{S}\right|}{\mathsf{\Delta }\mathrm{f}}_{\mathrm{l}}}$$

(44)

$${\mathsf{\Delta }\mathrm{f}}_{\mathrm{m}}=\left|\mathrm{f}\left({\mathrm{u}}_{\mathrm{m},\mathrm{g}}\right)-\mathrm{f}\left({\mathrm{x}}_{\mathrm{m},\mathrm{g}}\right)\right|$$

(45)

The amount of fitness improvement ∆fmis used in order to influence the parameter adaptation (S refers to either SCR or SF). As MCR is updated, if MCR,k,g = -1 or max(SCR) = 0 (i.e., all elements of SCR are 0), MCR,k,g+1 is set to -1. Thus, if MCR is assigned the terminal value -1, then MCR will remain fixed at -1 until the end of the search. This has the effect of locking CRi 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 has been shown to work well in conjunction with PSR methods. For the development of the present optimization algorithm, it has been decided to incorporate the continuous adaptive population reduction (CAPR) method. The CAPR method gradually reduces the population size according to the change of gradient of the fitness value [25].

$${\mathrm{N}\mathrm{P}}_{\mathrm{g}+1}=\left\{\begin{array}{cc}\sqrt[\mathsf{\gamma }]{{\mathsf{\Delta }}_{\mathrm{g}}/{\mathsf{\Delta }}_{\mathrm{g}-1}}& 0<{\mathsf{\Delta }}_{\mathrm{g}}/{\mathsf{\Delta }}_{\mathrm{g}-1}<1\\ {\mathrm{N}\mathrm{P}}_{\mathrm{g}}& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.$$

(46)

$${\mathrm{N}\mathrm{P}}_{\mathrm{g}+1}=\left\{\begin{array}{cc}{\mathrm{N}\mathrm{P}}_{\mathrm{g}+1}& {\mathrm{N}\mathrm{P}}_{\mathrm{g}+1}>{\mathrm{N}\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\\ {\mathrm{N}\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}& {\mathrm{N}\mathrm{P}}_{\mathrm{g}+1}\le {\mathrm{N}\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\end{array}\right.$$

(47)

where

$${\mathsf{\Delta }}_{\mathrm{g}}=\frac{{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{g}}\left({\mathrm{x}}_{\mathrm{g}}\right)-{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{g}}\left({\mathrm{x}}_{\mathrm{g}-1}\right)}{{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{g}}\left({\mathrm{x}}_{\mathrm{g}}\right)},{\mathsf{\Delta }}_{\mathrm{g}-1}=\frac{{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{g}}\left({\mathrm{x}}_{\mathrm{g}-1}\right)-{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{g}}\left({\mathrm{x}}_{\mathrm{g}-2}\right)}{{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{g}}\left({\mathrm{x}}_{\mathrm{g}-1}\right)}$$

(48)

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

Another criterion to consider when wanting to increase the performance of algorithms based on Differential Evolution is the way to generate 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 among the most popular sampling techniques in computer experiments thanks to its simplicity and projection properties with high-dimensional problems. 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 [26].

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

• Exhaustion-based criteria: Due to limited computational resources optimization run might be terminated after a certain generation, number of objective function evaluations or CPU time. Commonly, 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 is not able to stop the run.
• Distribution-based criteria: For DE algorithms all individuals converge to the optimum eventually. Therefore, it can be concluded that convergence is reached when the individuals are close to each other. Because is assumed that the optimum is not known as for the reference criterion, the distances between the population members are examined. This type of criterion can be applied in the design space or in the 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 of the ways to speed up the calculation process is by using parallel computing strategies. The strategy that was proposed to be used is Lightweight Pipelining (LP) [28]. The pipelining process helps in providing an easy approach in downloading and using the models on-demand. It helps in parallelization which means different jobs can be run parallelly also it reduces redundancy and helps to inspect and debug the data flow in the model. Some of the features that pipelines provide are on-demand computing, tracking of data and computation, inspecting the data flow, etc.

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. The coding of this algorithm was performed based on the Python 3 language mainly, in a GNU/LINUX environment.

Algorithm 4: OpenVINT algorithm
// Initialization phase
Input the design constants d, V∞, ρ, ν, a, Tmin, $\overline{r}$, NS;
Input the design variables intervals [Xmin, Xmax];
Input the optimization conditions G, NP, NPmin, ε, U*, γ, H, p, g=1;
Initialize of metrics;
Initialize population Pg with LHS;
//Parallelized loop by joblib
for i = 1 to NP do
	Apply Algorithm 1 for xi,g;
//Parallelized loop by joblib
for i = 1 to NP do
	Get T(xi,g), W(xi,g), ηd(xi,g), ηs(xi,g) with Algorithm 2;
Get ψ(xg) with (2) and L(xg) with (1);
Update U*;
Save Lavg(xg);
Save data of generation g;
Set all values in MCR, MF to 0.5;
Archive A = ∅;
k = 1;
// Main loop
for g = 1 to G do
	SCR = ∅, SF = ∅;
	for i = 1 to NP do
		ri = select from [1, H] randomly;
		Get CRi,g;
		Get Fi,g;
		Get mutation vector vi,g;
		Get trial vector ui,g;
	//Parallelized loop by joblib
	for i = 1 to NP do
		Apply Algorithm 1 with ui,g;
	//Parallelized loop by joblib
	for i = 1 to NP do
		Get T(ui,g), W(ui,g), ηd(ui,g), ηs(ui,g) with Algorithm 2;
	Get ψ(ug) with (2) and L(ug) with (1);
	Update U*;
	for i = 1 to NP do
		if L(ui,g) ≤L(xi,g) then
			xi,g+1 = ui,g;
			xi,g→ A;
			CRi,g→ SCR, Fi,g→ SF;
		else
			xi,g+1 = xi,g;
	Update memories MCR and MF with Algorithm 3;
	Save Lavg(xg+1);
	if g≥ 3 then
		Get Δg and Δg-1 with (48);
		Get NPg+1 with (46);
		if NPg+1<NPmin then
			Apply (47);
		(NPg – NPg+1)-th worst vectors → A;
		Remove the (NPg – NPg+1)-th worst vectors from Pg+1;
	Save data of generation g+1;
	if |Lavg(xg+1) – Lopt(xg+1)| ≤ε then
		break;
	k++;
Print metrics plots;
Output xopt, L(xopt);
Drawing the optimal propeller in point clouds;

3. Study Case

To evaluate the performance of the OpenVINT algorithm, a test was conducted to obtain the optimal design of a propeller used for an engine of a fixed-wing aircraft.

The flow characteristics for this study are V_∞, 25 m/s; ρ, 1.225 kg/m³; ν, 0.000014607 m²/s; and a, 340.294 m/s. A minimum permissible thrust of 7.5 N was considered for the optimization process. The intervals of the design variables are shown in Table 1.

For this optimization case the following parameters of the optimization algorithm will be taken:

• in real optimization problems it is considered to use at least 50 individuals in the initial population, and as a minimum population 10 individuals were considered;
• the stop conditions contemplated for this case were, a maximum number of evaluated generations of 200, an ε value of 1 W to fulfill the condition indicated in line 48 of algorithm 4;
• a γ factor of 50 was used in equation 64;
• finally, an initial U*-value of 100 W was considered.

4. Results

To find the optimal propeller geometry, the algorithm needed to evaluate 127 generations to meet one of the stopping conditions. The penalty function was evaluated 5742 times. Figure 6 shows the evolution of the optimization process.

The optimal propeller can generate a thrust of 7.505 N, requiring a power of 226.08 W, with a dynamic efficiency of 83%. The optimal vector is shown in Table 2. The next figure shows the geometric characteristics of the propeller.

Geometric characteristics of the optimal propeller are shown in Figure 7.

The general view of the propeller and the shape of the airfoils in each section in $\overline{r}$ of the propeller blade are shown in Figure 8.

4. Discussion

Good agreement was shown by comparing the calculations of the aerodynamic characteristics of airfoils in Xfoil with the results of solving the Navier-Stokes equations using the control volume method and experimental data. This allows fast discrete eddy Xfoil calculations with viscosity and compressibility corrections to be used for propeller design calculations.

The cross-sectional calculation method implemented for propeller design showed good agreement with experiment and CFD modeling in terms of thrust and torque (Figure 5a,b) and provided a quick time for assessing aerodynamic characteristics.

The optimization of propeller parameters was based on the differential evolution method. It is shown that 120 generations of 60 individuals are sufficient to determine 27 design variables.

The solution to the demonstration problem of designing an aircraft propeller with a required thrust of 7.505 N with a diameter of 300 mm and a free-stream speed of 25 m/s yielded a combination of design parameters that provides a propeller efficiency of 83% and requires 226 W of power at a rotation speed of 6896.5 rev/min.

Approbation of the proposed methodology in the design of real propellers of unmanned aerial vehicles has shown its effectiveness. The propellers obtained because of optimization have consistently high efficiency. This methodology has high performance and does not require large computing power.

5. Conclusions

This paper describes the developed methodology for designing air propellers. The novelty of the technique is the selection of the aerodynamic profile for each cross-section of the propeller blade by using Bezier curves for accuracy, flexible and robust setting of the geometry of the blades. It is shown that to accurately describe a propeller, 27 parameters are sufficient, of which 16 are for describing the profile, 10 are for the shape of the blade, and 1 is for specifying the kinematics (rotation speed of the propeller).

The use of discrete vortex and isolated section methods makes it possible to calculate the thrust and required propeller power in less than one second. Combined with the differential evolution method, this makes it possible to solve the optimization problem in less than 2 hours on a workstation with four Intel Core i5 cores.

The OpenVINT program has been developed in Python, which implements the presented methodology with the possibility of parallel calculations (standardly with four threads). To evaluate the performance of the OpenVINT algorithm, a test was conducted to determine the optimal design of the propeller used for the engine of an unmanned aerial vehicle. The algorithm makes it possible to obtain propellers with parameters that provide an efficiency of up to 80% under given requirements, which is a good indicator.

The program that implements the technique also allows you to automatically obtain a geometric description of the obtained result in the form of point clouds, from which it is possible to quickly construct a three-dimensional geometric model in any of the existing CAD systems (see Figure 8) for the purpose of its further use in subsequent calculations in other programs or for the development of equipment for the manufacture of propellers using any available technologies.