4.1 Modelling overview

BEMT uses both blade element theory and one-dimensional momentum theory to calculate through an iterative process the induced velocities (represented by the induction factors) on each station of a discretised blade. Blade element theory deals with relative velocities and angles, while momentum theory deals with the relationship between the velocities at different points of the flow. A detailed description of both methods is presented by Glauert in Ref.^{(Reference Glauert23)}.

Figure 2.

Definition of blade element.

Figure 2 shows a blade element, and graphically explains the concept of flow angle, $\phi$ , which is the angle between the plane of rotation and the local relative velocity, which is dependent on the axial airspeed and the local tangential speed, the latter being variable along the blade radius. Induction factors can be defined as follows: the axial induction factor $a_0$ is the relative increment of axial speed from upstream to the plane of rotation, while the tangential induction factor $a_1$ is the relative decrement of the local tangential speed on each blade element. These changes to the axial and tangential velocities are associated with induced velocities and affect the flow angle of each blade section. According to Fig. 2, the flow angle can then be expressed as shown in Equation (10). The induced velocity is assumed to be parallel and opposite to the local lift vector, which is generally true since the aerofoil operates in its linear region; otherwise, this assumption is not strictly correct but still provides enough accuracy. According to Fig. 2, the flow angle can then be expressed as shown in Equation (10).

(10) \begin{equation} \phi = \tan ^{-1} \left[ \frac{V_\infty(1+a_0)}{\Omega r (1-a_1)} \right] = \tan ^{-1} \left[ \frac{\Omega r a_1}{V_\infty a_0}\right] \end{equation}

The method divides the blade into discrete independent elements with an associated aerofoil whose lift and drag coefficients are known for a range of angles of attack. Knowing the flow angle of each blade element and the local pitch, $\beta$ , defined as the angle between the local chord and plane of rotation as shown in Fig. 2, the angle-of-attack of each blade element can be estimated through Equation (11), then the corresponding lift and drag coefficients can be computed from a known database; the lift and drag coefficients of each section are decomposed into normal and tangential components as shown in Equations (12) and (13).

(11) \begin{equation} \alpha = \beta - \phi \end{equation}

(12) \begin{equation} c_t = c_l \cos \phi - c_d \sin \phi \end{equation}

(13) \begin{equation} c_q = c_l \sin \phi + c_d \cos \phi \end{equation}

From the calculated normal (thrust) and tangential coefficients, the normal and tangential force contribution can be computed for each blade element, forming a force distribution which when integrated allows an estimation of the total thrust and torque, as shown in Equations (14) and (15).

(14) \begin{equation} T=B \int_{R_{hub}}^{R} \frac{1}{2} \rho V_{rel}^{2} c_t c \text{d}r \end{equation}

(15) \begin{equation} Q=B \int_{R_{hub}}^{R} \frac{1}{2} \rho V_{rel}^{2} c_q c r \text{d}r \end{equation}

From momentum theory, the blade is analysed as an actuator disk divided into discrete rings, each representing a blade element. This analysis provides a different way to estimate the normal and tangential force contributions, which when combined with that calculated by blade element theory, gives expressions to calculate both the normal and tangential induction factors, as shown in Equations (16) and (17).

(16) \begin{equation} a_0 = \left( \frac{4F \sin ^2 \phi }{\sigma c_t} -1 \right)^{-1} \end{equation}

(17) \begin{equation} a_1 = \left( \frac{4F \sin \phi \cos \phi }{\sigma c_q} +1 \right)^{-1} \end{equation}

where $\sigma$ is defined as the solidity factor (see Equation (18)), the ratio between the actual solid area of the blade element and the area swept out by that element.

(18) \begin{equation} \sigma = \frac{B c}{2\pi r} \end{equation}

AN iterative process can then be carried out for each blade section by assuming both induction factors to be zero as an initial guess, then the flow angle can be calculated through Equation (10). Next, the angle-of-attack is estimated as defined in Equation (11), for which the lift and drag coefficients can be computed from the known blade section aerofoil. The lift and drag coefficients are transformed to normal and tangential coefficients using Equations (12) and (13), which provides a new estimation of the induction factors through Equations (16) and (17). The process is then restarted with the new estimates of the induction factors; when their values converge, the normal and tangential force distributions can be computed and finally integrated to yield the thrust and torque. This method is implemented computationally to provide a propeller performance prediction tool for the optimisation algorithm. A flowchart for this computational BEMT tool is shown in Fig. 3.

Figure 3.

BEMT code flowchart.

4.2 Corrections

4.2.1 Drag coefficient correction

Lift and drag polars are obtained from CFD simulations and extrapolated through the Montgomerie method, but the CFD simulations should be carried out at specific conditions, that is, specific Reynolds number; the selection of the Reynolds number for the aerofoil simulations depends on the radial position of the aerofoil. Ideally, the aerofoil should be simulated at the Reynolds number correspond to its radial position and the operational condition of the propeller. Nevertheless, to keep the computational cost of this method low, the aerofoil is simulated at one specific Reynolds number, then the drag coefficient is corrected to account for the actual Reynolds number^{(Reference HernÁndez and Crespo24)}. The correction is defined in Equation (19):

(19) \begin{equation} c_{d_{\mathrm{Re}}} = c_d \left( \frac{\mathrm{Re}_{ref}}{\mathrm{Re}} \right)^{0.2} \end{equation}

The lift coefficient is not corrected as the linear part of the lift curve does not change considerably with Reynolds number. The stall region is affected but a Reynolds correction for this has not been developed so far.

4.2.2 Prandtl’s tip and hub-loss approach

The fact that there are a finite number of blades of finite length affects the force distribution, modifying it due to three-dimensional and tip effects, as is the case for the finite wings of an aircraft. To account for this phenomena, Prandtl’s tip loss approach is implemented. There is a Prandtl’s tip-loss factor associated with each blade station, which modifies the local loads, namely the normal and tangential forces; this factor varies along the radius and is close to unity near the root but decreases to zero towards the tip, meaning that tip losses barely affect the loads near the root while at the tip they completely vanish. An approximate method to calculate Prandtl’s tip-loss factor for each section of a rotor is presented in Ref.^{(Reference Sriti and Khchine25)} and shown in Equation (20); this method developed for wind turbines is also valid for propellers.

(20) \begin{equation} F_{tip} = \frac{2}{\pi}\cos ^{-1} \left[ exp\left( -\frac{B(R-r)}{2r\sin \phi} \right) \right] \end{equation}

The wake shed by the hub also affects the blade loading distribution towards the root^{(Reference Rwigema26)}; the hub-loss correction factor also presented in Ref.^{(Reference Sriti and Khchine25)} is very similar to the tip-loss correction factor and is defined in Equation (21).

(21) \begin{equation} F_{hub} = \frac{2}{\pi}\cos ^{-1} \left[ exp\left( -\frac{B(r-R_{hub})}{2r\sin \phi} \right) \right] \end{equation}

Finally, to account for both the tip and hub losses, both correction factors are multiplied (see Equation (22)) to obtain the final correction factor to be applied at each blade station.

(22) \begin{equation} F = F_{hub}F_{tip} \end{equation}

4.2.3 Correction for rotational effects

The lift and drag coefficients of each blade section, which finally feed the BEMT method, are generally available as two-dimensional data, estimated through wind-tunnel tests, CFD simulations or software such as XFOIL, which assumes that the aerofoil is subjected to 2D flow; for the case of rotors, such as propellers or wind turbines, each blade section is subjected to rotational movement; this fact makes two-dimensional data inappropriate for rotor calculations because of phenomena such as centrifugal and Coriolis forces to which the boundary layer is subjected^{(Reference Hansen, SØrensen, Voutsinas and Madsen27)}. These rotational effects are associated with the three-dimensionality of the problem; one of the simple rotational corrections shown in Ref.^{(Reference Bak, Johansen and Andersen28)} is used in the current method, being a semi-empirical correction to the lift coefficient as a function of the station’s radial location, its chord and its local relative velocity, as shown in Equation (23). According to Ref.^{(Reference Bak, Johansen and Andersen28)}, the value of the coefficient b is $3.1$ for rotors; but for this study, the factor takes a value of $1.5$ , as selected by McCrink in Ref.^{(Reference Mccrink and Gregory21)}, p. 3 in more similar work. Note that the term $c_{l_{pot}}-c_{l_{2D}}$ accounts for the stall and viscous effects on the correction.

(23) \begin{equation} c_{l_{3D}} = c_{l_{2D}} +b\left( \frac{c}{r} \right)^2(c_{l_{pot}}-c_{l_{2D}})\left( \frac{\Omega r}{V_{rel}} \right)^2 \end{equation}

4.2.4 Mach correction

The lift and drag coefficients are affected when the compressibility of the fluid becomes more significant, which is when the relative velocity of a given blade element reaches a Mach number higher than $0.3$ ; this has special importance for the elements towards the tip of the blade when it is operating at high rotational speed. That is why the lift and drag coefficients feeding the BEMT method are corrected to account for compressibility effects through the Karman–Tsien correction^{(Reference Tsien29,Reference von Karman30)} . The lift coefficient correction is shown in Equation (24), which is also valid for the drag coefficient.

(24) \begin{equation} c_{l_{M}} = \frac{c_l}{\sqrt{1-M^2}+\frac{\displaystyle c_l M^2}{\displaystyle 2 + 2\sqrt{ 1-M^2}}} \end{equation}

5.1 Motivations

Commonly, numerical potential theory via panel methods, such as XFOIL (which also uses a boundary-layer model), are implemented in the BEMT algorithm to predict the aerodynamic coefficients^{(Reference MacNeill, Verstraete and Gong19,Reference Mccrink and Gregory21,Reference Gu, Cai, Zhou, Li, Shen and Zhang31,Reference Wall32)} . Although this results in good approaches with low computational cost, low Reynolds conditions and the post-stall regime are not well predicted by such codes, furthermore, XFOIL overpredicts the maximum lift coefficient and underpredicts the minimum drag^{(Reference Coder and Maughmer33)}.

The lift and drag coefficients in the post-stall regime can be predicted by extrapolation techniques (Viterna and Janetzke^{(Reference Viterna and Janetzke34)}, Montgomerie^{(Reference Montgomerie13)}, Lindenburg^{(Reference Lindenburg35)} and Spera^{(Reference Spera36)}), but the angle at which the aerofoil achieves the maximum lift is not predicted by these extrapolations, thus the maximum lift coefficient inaccuracy remains if such methods are used to extrapolate XFOIL predictions.

On the other hand, CFD methods have been demonstrated to provide more reliable results than XFOIL^{(Reference Barrett and Ning37-Reference Morgado, Vizinho, Silvestre and PÁscoa39)}, especially in the non-linear regime and low Reynolds conditions, as long as a proper turbulence model and mesh are employed. Considering that correct BEMT predictions decisively depend on the aerofoil performance data^{(Reference Koch40)}, CFD methods are adopted here to predict the aerofoil linear and stall regions. Then, to keep the computational cost relatively low, the post-stall region is approximated through the Montgomerie extrapolation method^{(Reference Montgomerie13)}, which models the aerofoil post-stall region as a flat plate, giving acceptable results since the maximum lift coefficient is well estimated. This method presents a smoothed transition model between the stall and post-stall regions, making it the most reliable of the extrapolation methods^{(Reference Bianchini, Balduzzi, Rainbird, Peiro, Graham, Ferrara and Ferrari41)}. This CFD–Montgomerie combination offers accurate pre-stall and stall aerofoil performance prediction thanks to the robustness of CFD methods and quite good estimation of the post-stall region, which is not extremely decisive for an accurate propeller performance prediction.

5.2 Turbulence model selection

There are three categories of models that can be implemented to solve a turbulent flow: Large-Eddy Simulation (LES), Direct Numerical Simulation (DNS) and RANS^{(Reference AraÚjo, Silva, Santos, Corte Real, Menezes, Rocha, Rodrigues, Borges, Costa and Bezerra42)}. The former two are not usually employed due to their high costly, even more so when the aim is to implement it in several simulations as in the current case. Among the wide range of RANS models available, the transition Shear Stress Transport (SST) model is selected. This four-equation turbulence model is based on the $k-\omega$ SST model^{(Reference Menter43)}, but additionally, two other equations are added: one for the intermittency, $\gamma$ , and another for the laminar–turbulent transition with $\mathrm{Re}_{\theta t}$ criteria, which links empirical transition data with the intermittency equation^{(Reference Wang, Ingham, Ma, Pourkashanian and Tao44)}. This four-equation model is also termed the $\gamma- \mathrm{Re}_{\theta} -$ SST model^{(Reference Menter, Langtry and Völker45)}, being commonly used as it provides reliable results for the linear and stall regime, including low Reynolds number condition^{(Reference Munduate and Ferrer38,Reference AraÚjo, Silva, Santos, Corte Real, Menezes, Rocha, Rodrigues, Borges, Costa and Bezerra42,Reference GÜnel, KoÇ and Yavuz46)} .

5.3 Grid convergence and simulation settings

Due to the large number of simulations required to cover the linear and stall regimes for each selected aerofoil, it is important to keep the computational costs low, hence a grid convergence study was carried out using a coarse, medium and fine mesh.

The domain is a structured circleFootnote 2 with a radius of 100 chords. This kind of O-ring type of grid allows the use of only one mesh for all angles of attack, imposed by specifying the flow direction at the outer inflow boundary of the grid, as also done in Refs^{(Reference Morgado, Vizinho, Silvestre and PÁscoa39,Reference GÜnel, KoÇ and Yavuz46)} . The aerofoil surface has a no-slip boundary condition. The distance from the aerofoil to the first layer is calculated by flat-plate boundary-layer theory^{(Reference White47)}, p. 467 as $8.524 \cdot 10^{-5}$ m for the three meshes with the aim of keeping $y^+ \leq 1$ at the first grid point off the surface, which is recommended for better CFD turbulent predictions^{(Reference Salim and Cheah48)}.

The aerofoil tested is a NACA4412 at Re=250,000. The simulations are performed using the ANSYS FLUENT steady-state solver. A density-based solver is applied. The external circumference is defined as the far-field pressure with M = 0.01074; this kind of boundary requires an ideal gas law to be implemented. A Second-Order Upwind (SOU) scheme is used for spatial discretisation of the flow, turbulence and momentum equations^{(Reference Morgado, Vizinho, Silvestre and PÁscoa39)}, p. 210. The conservative equations are solved by an implicit formulation with the Advection Upstream Splitting Method (AUSM) flux type, and finally, the pseudo-transient option is implemented with a timescale factor of $0.001$ to achieve fast convergence.

The main reason why CFD methods are employed is that they represent an adequate approach for stall phenomena, thus the maximum lift coefficient is a proper criterion to compare the changes in predictions based on the grid resolution. The maximum lift coefficient for each mesh and the grid data is presented in Table 1.

Table 1

Grid independence data and results

According to the data presented in Table 1, the medium refinement was selected to run the simulations due to the small difference in its maximum lift coefficient prediction compared with the fine refinement. The detail of the medium refinement mesh is shown in Fig. 4.

5.4 CFD–Montgomerie validation

The NACA0012 aerofoil was selected to validate the aerodynamic predictions of the CFD–Montgomerie methodology, where the transition between the CFD and Montgomerie approaches occurs at an angle-of-attack of 10 $^{\circ}$ , which is the stall angle predicted by CFD. Results are shown in Fig. 5. This method predicts the 2D aerofoil data for further 3D correction.

Figure 4.

Medium mesh overview for NACA4412.

Figure 5.

CFD–Montgomerie validation for NACA0012 at $\mathrm{Re}=360,000.$

From the data shown in Fig. 5, the drag and lift coefficient estimations have a Mean Absolute Percentage Error (MAPE) of 22.98% and 16.55% respectively. This does not represent a critical problem since the pre-stall regime shows a lift MAPE of 4.43% and most of the aerofoils on the blade lie in this regime.

7.1 Number of blades

The argument about the optimal number of blades requires a brief explanation about the practical changes when the number of blades is increased. As Equation (8) shows, the propeller efficiency depends on a specific condition of the airspeed and propeller angular speed, so to enable a proper comparison, the same conditions must be considered, including the operating thrust.

Furthermore, the thrust developed as well as the power consumption are related to the blade’s area. Hence, an increase in the number of blades must be countered by a reduction of the diameter or chord length to maintain the same thrust, but a reduction of both the diameter and chord length causes a decrease in the Reynolds number, consequently affecting the performance of the section aerofoils, essentially due to a decrease of the boundary-layer thickness and the movement of the transition point^{(Reference Anderson54)}, pp. 375–377.

In Ref.^{(Reference Dai, Quan, Ren and Cai53)}, a model corrected with experimental data, which shows the relation between the efficiency and the number of blades, is developed:

(25) \begin{equation}\eta_{P}\propto \frac{1}{B^{(\alpha_m - \alpha_t)}}\end{equation}

where $\alpha_m$ , $\alpha_t $ are correction coefficients calculated from experimental data. As Equation (25) shows, a decrease in the number of blades causes an increase in the efficiency because $\alpha_t < \alpha_m $ , although the use of one blade is not recommended for balance reasons. Thus, the optimal number of blades is $B=2$ . Nevertheless, if there is a diameter constraint, it is possible that a three-bladed propeller may be the only way to fulfil the thrust requirements, thus the algorithm allows two- and three-bladed propellers.

7.2 Chord and pitch distribution

Prandtl and Betz developed a vortex theory for application to propellers^{(Reference Betz15)}, analogous to lifting line theory for a wing. Vortex theory models the wake shed by the propeller as a helical vortex sheet (Fig. 8) that induces velocities along the blades, thus changing the local angle-of-attack at each blade section^{(Reference Larrabee55)}.

Figure 8.

Propeller vortex model.

The theory states that there is a certain circulation distribution along the blade which minimises the loss of kinetic energy, being analogous to the elliptical lift distribution of a wing that minimises the induced drag. Just as for wings the condition for minimum induced drag is a constant induced angle along the span, for propellers the ideal condition is a constant horizontal velocity at each helical wake location. There are certain pitch and chord distributions that achieve this condition while keeping the other geometrical parameters fixed. The method for the estimation of such a distribution developed by Prandtl and Betz was improved by Goldstein^{(Reference Goldstein14)}, and this is selected to find the optimal pitch and chord distributions for the optimisation algorithm. This optimisation is performed for a specific diameter and a specific operational condition, namely airspeed, rotational speed and target thrust.

First, a local dimensionless velocity is defined as shown in Equation (26):

(26) \begin{equation} x = \frac{\Omega r}{V_\infty} \end{equation}

Then, a new variable associated with the local circulation, G, is defined and then expressed as a function of x as shown in Equation (27).

(27) \begin{equation} G = \frac{B \Omega \Gamma}{2\pi V_\infty V'} = \frac{F_{tip}x^2}{1+x^2} \end{equation}

For feasibility of modelling, a more convenient definition of the thrust coefficient is established in Equation (28).

(28) \begin{equation} C_t = \frac{2T}{\rho V^2_{\infty} \pi R^2} \end{equation}

This form of the thrust coefficient is used only for pitch and chord optimisation.

Then, $\zeta$ is defined as the ratio between the helical vortex horizontal speed and the airspeed. An expression for it was developed by Goldstein^{(Reference Goldstein14)} and is shown in Equation (29).

(29) \begin{equation} \zeta = \frac{V'}{V_{\infty}} =\frac{I_1}{2I_2}\left(1-\sqrt{1-\frac{4 I_2 C_t}{I_1^2}} \right) \end{equation}

with $I_1$ and $I_2$ defined as shown in Equations (30) and (31).

(30) \begin{equation} I_1 = 4 \int_{0}^{1}G\left(1-\frac{1}{x(L'/D')_{max}}\right)\xi d \xi \end{equation}

(31) \begin{equation} I_2 = 2 \int_{0}^{1}\frac{G}{x^2+1}\left(1-\frac{1}{x(L/D)_{max}}\right)\xi d \xi, \end{equation}

defining $\xi$ as the local dimensionless radius as shown in Equation (32).

(32) \begin{equation} \xi = \frac{r}{R} \end{equation}

Now, the optimal chord for each blade section can be calculated through Equation (33) for the corresponding values of G and x and a given local lift coefficient, being that corresponding to the angle-of-attack at which the maximum lift-to-drag ratio is achieved.

(33) \begin{equation} c_{opt} = \frac{4\pi R\zeta \lambda G }{c_{l_{opt}}\sqrt{1+x^2}} \end{equation}

with $\lambda$ being the dimensionless speed at the blade tip, as Equation (34) shows.

(34) \begin{equation} \lambda = \frac{V_\infty}{\Omega R} \end{equation}

Additionally, the local flow angle that minimises the induced losses is estimated as shown in Equation (35).

(35) \begin{equation} \phi_{opt} = \tan ^{-1}\left[\frac{\lambda}{\xi}\left(1+\frac{1}{2}\zeta \right) \right]\end{equation}

Now, knowing the flow angle and the angle-of-attack at which the local aerofoils present their best lift-to-drag ratio, the optimal pitch can be estimated through Equation (36).

(36) \begin{equation} \beta_{opt}=\phi_{opt} + \alpha_{opt} \end{equation}

7.3 Aerofoil selection

The optimisation algorithm evaluates a list of aerofoils along the propeller, which are selected by taking into account the kind of aerofoil employed for each section of a small propeller^{(Reference Tjhai52,Reference Reuben, Hiraki and Shohei56,57)} . The aerofoil varies along the propeller. At the root there is a low Reynolds aerofoil with a strong camber, while at the tip a less cambered aerofoil is employed^{(Reference Glauert23)}. This decrement in camber along the radius is implemented due to structural criteria. Once the camber filter is applied, the criterion to select the aerofoils is the best lift-to-drag ratio or aerodynamic efficiency estimated from Ref.⁽⁵⁷⁾. Tree different aerofoils for three different sections of the blade (root, mid and tip) are selected for the optimisation process, and the maximum camber and thickness for each aerofoil are presented in Table 3.

The list presented in Table 3 is aimed at small-scale propellers and their relatively low Reynolds number regime. In general, the Reynolds numbers expected for the intended propeller design dictate the preliminary aerofoil selection. The Reynolds number at which the aerofoil data will be estimated for the tip, middle and root aerofoil simulations is selected to be 530,000, 130,000 and 90,000 respectively.

Table 3

List of selected aerofoils

7.4 Algorithm

The algorithm includes three main subroutines: the Goldstein optimisation, the BEMT algorithm and the electric motor model. The first step in the optimisation algorithm is to set the design parameters and the range for the variables to optimise, which are both inputs for the Goldstein and BEMT methods: the cruise speed, a certain flight altitude that allows the calculation of air density, sound speed and viscosity through the standard atmosphere and Sutherland’s law, the target thrust or $D_{wake}$ from Equation (1), the number of blade sections and the ranges for the design variables such as d, B and $\Omega_{guess}$ , which is a guess value for the vortex algorithm due to its prediction discrepancy compared with BEMT results, this value being very similar to the operational rotational speed. Finally, the root, mid and tip aerofoils are selected from the list shown in Table 3, and the position of the middle aerofoil is also varied for the optimisation. Additionally, the lift and drag coefficients are interpolated between the known aerofoil performance estimated from the CFD–Montgomerie method to account for the aerofoil shape transition.

The next step is the Goldstein optimisation, whose outputs are the optimal chord and pitch distributions for the design conditions, including $\Omega_{guess}$ . Then, the BEMT algorithm (Fig. 3) employs this data to calculate the propeller performance, whose useful outputs are the thrust and torque as functions of the rotational speed and $\eta_{P}(J)$ , then the rotational speed at which the target thrust is achieved is found. Subsequently, the operational advance ratio is calculated taking into account the estimated rotational speed to finally obtain the propulsive efficiency.

Finally, the electric inputs $i_o$ , $K_V$ , $R_a$ and $R_{ESC}$ are employed for the electric motor model. The torque is calculated from Equation (8), then the current consumption is found from Equation (4) and the voltage from Equation (6). The electric performance is summarised by i, v and $\eta_{E}$ calculated from Equation (7), and $\eta_{T}$ from Equation (9).

Some constraints are applied to discard solutions that exceed the maximum allowable power and voltage of the motor. Additionally, there is another internal constraint to restrict the propeller tip speed to Mach 0.7 (as used in Ref.^{(Reference Gur and Rosen58)}, p. 1344), with the aim of ensuring that the speed of sound is not exceeded at any point over the aerofoil. The solution with the minimum electrical power consumption, $P_{E}$ , is selected as the optimal design point to maximise the endurance. The complete algorithm is subjected to MATLAB’s particle swarm optimisation routine to find the optimum values of diameter, rotational speed, distribution of aerofoils and number of blades.

8.1 Optimisation case

The optimisation results are validated through a comparison with an existing optimisation case presented in Ref.^{(Reference Gur and Rosen59)} at sea level, whose main data and constraints are presented in Table 4.

Table 4

Data and constraints for the optimisation comparison

In addition to the data presented in Table 4, the algorithm is tested for a two- and three-bladed propeller with a hub radius of 10% of the total blade radius and a no-load current of 0A. The range to test the diameter is from 0.1m to 0.3m, $\Omega_{guess}$ from the value corresponding to J = 1 to the one correspondent to 0.7M at the blade tip, and finally the position for the middle aerofoil is tested from 40% to 70% of the radius length with three steps. The results of the optimal propeller parameters and the electric performance are presented in Table 5.

The optimal propeller designed in Ref.^{(Reference Gur and Rosen59)} theoretically achieved a power consumption of $P_E=56$ W. The current design predicts a minimum energy consumption of $P_E=54.2$ W with the same value for d and different values for variables as $\Omega$ , and chord and pitch distributions shown in Fig. 9(a) and (b), respectively. These distributions computed from the Goldstein algorithm are fitted using a third-degree polynomial to avoid too abrupt changes.

Table 5

Study case optimal propeller performance

Figure 9.

Propeller optimal chord and pitch distribution case.

Figure 10.

Impact of d and $\Omega$ on $P_E$ .

Based on the cases generated by the algorithm, some important trends related to power consumption can be studied. The main tendency reveals that, for a given target thrust, an increase in the propeller diameter and operational rotational speed results in lower power consumption due to the consequent increment of tangential speed, resulting in a higher operational Reynolds number, a condition taht improves the aerodynamic efficiency of section aerofoils. This tendency is shown in Fig. 10 for a two-bladed propeller with the aerofoils selected for the optimal case described in Table 5.

Additionally, for different generated cases, the behaviour of the propulsive, electric and total operational efficiency is analysed for a specific range of rotational speeds, as shown in Fig. 11. The data show how the optimal operational condition for the propeller (the point with the highest $\eta_P$ ) is not close to the optimal operational point of $\eta_E$ . The trend shows that the motor achieves better values of electric efficiency at high values of rotational speed. On the other hand, the peak efficiency of the propeller is achieved at a different value of rotational speed, which is the main reason to use a gearbox for this kind of system. Therefore, to improve the performance of the propulsion system, it is not enough to optimise either the propeller or the motor operational condition alone; rather, it is important to find an equilibrium point where both the propeller and electric system operate at an acceptable performance level while keeping the overall performance as high as possible.

Figure 11.

$\eta_P$ , $\eta_E$ and $\eta_T$ for different design rotational speeds ( $d=0.262$ m, $B=2$ ).

From the study case, a tendency chart is constructed to show how the number of blades affects the operational propulsive efficiency at different operational rotational speeds, as shown in Fig. 12. It is shown that, for high rotational speeds, higher efficiencies can be reached because of the higher operational Reynolds number of the whole blade and its positive effect on aerofoil performance; additionally, the tendency shows that, for fixed diameter, at lower rotational speeds, the advantage of using a three-bladed propeller is more marked than at higher rotational speeds, at which even a two-bladed propeller appears to offer slightly higher performance.

Figure 12.

$\eta_P$ at different operational rotational speed for 2 and 3 blades ( $d=0.17$ m).

8.2 CFD comparison

To evaluate the performance of the optimal propeller designed in the case study presented above, numerical analyses are performed using Autodesk CFD software, which employs the finite element method for the discretisation. The aim is to compare the reliability and validity of the BEMT method with the CFD results. The simulation employs a rotatory region which is a cylindrically shaped region around the propeller. The external boundary has a velocity inlet boundary of 10m/s and an outlet face with 0Pa gauge pressure. A structured quadrilateral mesh with 15 layers is employed along the blade surface with a growth rate of $1.2$ . One blade of the propeller and the CFD flow interaction are shown in Fig. 13.

Figure 13.

Optimal propeller CFD study.

The SST k– $\omega$ model is selected for the turbulent flow with a time step of $0.00005$ s and 20 inner iterations. The result shows a total force of $T=3.28$ N for two blades, without considering the hub drag, which is not predicted by the BEMT method. The prediction presents a relative error of $6.28$ %.