Advanced Scale-Propeller Design Using a MATLAB Optimization Code

1. Introduction

1.1. Rationale

A propeller is a propulsion device that exerts a force on a large mass of air. The earliest versions of the propeller have been around for centuries, with Archimedes developing a rotating screw system, which could be used to draw water from deep wells. Since their first use in the 19th century, hydrodynamic propellers have been the propulsion system of choice for countless marine vehicles. However, it wasn’t until the 20th century when the Wright brothers used the first aerodynamic propeller as the propulsion system for their historic aircraft, the Wright Flyer. The Wright brothers were forced to develop their own propeller design method, due to the lack of existing aerodynamic propeller theory available to them. The method they used is now known as blade element theory, and it is still used to design propellers today. In simple terms, blade element theory assumes that the airflow over each section of a propeller blade is two-dimensional. Whilst in reality there will be aerodynamic interference between blade sections, it is assumed to be negligible. After several design iterations, the Wright brothers achieved a maximum propulsive efficiency of 78% by 1910. In comparison, modern large-scale propellers can now achieve maximum efficiencies of up to 90% [1].

However, this very high efficiency figure is infeasible for small scale propellers, such as those used on small FPV drones. Data analysis of the UIUC Propeller Database 1 & 4, by the authors, as shown in Figure 1, has shown that efficiencies of these small propellers (7-21” diameter) is typically between 45-75% [2], with a lot of variability between manufacturers, even with the same diameter and pitch. The reasons for this efficiency loss will be discussed in detail later in this paper, but relate to low Reynolds numbers, choice of airfoil and physical geometry. Research into this topic is limited to a few studies [2,3,4], meaning that there is ample scope for development in this important area of research.

1.2. Aim of the Study

The main aim of this work is to analyse the efficiency of scale propellers as used on small drones (see Figure 2) and to develop a MATLAB program which can accurately and easily optimise propeller geometry at operating conditions specified by the end user to improve efficiency and thereby increase endurance.

1.3. Background

High propulsive efficiencies can be achieved on propellers used on full-scale crewed aircraft, such as the C-130J Super Hercules, and there are two main reasons why this is the case. Firstly, the propellers tend to have complex variable pitch mechanisms, which allow the propeller to operate at peak efficiency under different flight conditions. Such mechanisms are costly, complex and too heavy to include within small propellers for attritable uncrewed aircraft. Variable pitch mechanisms are also required to control helicopter rotors, which is why these rotors tend to have a constant angle of twist along the blade.

Secondly, full-scale aircraft propellers tend to have large diameters, which leads to greater chord lengths. As will be explained later in this study, large diameter propellers operate in higher Reynolds (Re) number regimes, which enable greater efficiencies to be achieved, in contrast to smaller propellers which operate at lower Re numbers. Therefore, creating propeller geometries that permit higher efficiencies to be reached, whilst not inducing large mass and cost penalties, would lead to significant improvements in performance.

A previous research project at the University of Southampton resulted in a MATLAB code being developed which optimised propeller geometry for given operating conditions [4]. This work was used as the basis for the Mejzlik propeller online performance program [5]. Another research project examined the low Reynolds number propeller problem, but instead focused on using overlapping propellers as a solution [6]. However, neither of these projects attempted to design a single propeller which could operate at high efficiencies in low Reynolds number regimes.

1.4. Literature Review

The earliest analytical work on propellers was performed by Rankine and Froude in the 19th century. The theory developed was intended for use on marine propellers. The theory they developed later became known as Actuator Disk Theory (ADT)—or axial momentum theory. Whilst this doesn’t include any detailed propeller aerodynamics, this theory allows designers to develop rough performance and sizing requirements.

• ADT is based on a number of assumptions:
• The rotor is modelled as an actuator disk, which adds momentum and energy to the air;
• There is no inflow or outflow through the wake boundary;
• The flow is steady, incompressible and inviscid;
• The flow is one dimensional and uniform;
• The flow is uniform through the disk and wake;
• The disc does not impart any swirl to the flow.

Using inviscid flow theory and taking a force balance across the actuator disk, the thrust can be easily calculated. For a given thrust, the induced axial velocity can be calculated according to Eqn. 1.

$$w_{axial}=\sqrt{\frac{V^2}{4}+\frac{T}{2\mathsf{\rho }A}}-\frac{V}{2}$$

(1)

1.5. Historical Perspective

In the minutes of the (1890) Proceedings of the Institution of Civil Engineers [7], the use of propellers for marine applications was discussed, and it soon became clear that the theory they were developing was applicable to aeronautics as well.

Drzewiecki [1892] is credited with the development of Blade Element Theory (BET) [8]. BET is a method used to analyse propeller performance that, unlike actuator disk theory, considers the propeller geometry and aerodynamic properties. In blade element theory it is assumed that the flow over a blade section is two-dimensional. This means that the flow over a given blade section is unaffected by the flow over neighbouring sections—in this analysis, the effective velocity is simply the result of the axial velocity and the rotational velocity. The propeller blade is broken down into multiple sections, and the forces acting on each of these blade sections is determined. By integrating over the length of the blade, the overall blade forces can be developed.

Weick (1925) in his NASA technical note entitled: Simplified Propeller Design for Low-Powered Airplanes [9] detailed an empirical design method for propellers. The design process is accomplished through the use of charts and basic formulae and involves no detailed calculations. The data required for this design method is the aircraft airspeed, the brake horsepower and the revolutions per minute of the engine, and the density of air. From this data, the performance coefficient, advance ratio and efficiency of the propeller are found empirically. The article also educates on what materials are capable of withstanding the centripetal forces at different propeller diameters and rotational speeds.

In his 1935 paper, entitled Airplane Propellers [10], Glauert formulated a method for analysing arbitrary propeller designs. The method is a combination of axial momentum theory and blade element theory, and includes corrections for momentum loss due to radial flow. The method is limited, as the contraction of the propeller wake is not considered. As a result, the theory is accurate for low disc loading (small thrust per unit disc area), but is only adequate for estimating performance at high disc loadings.

Propeller Analysis from Experimental Data [11] is a 1940 report (Stickle & Crigler) in which propeller performance is analysed. The data used in this analysis was obtained in a 20-foot wind tunnel with a 4-foot diameter propeller, but the results apply to all propellers. Firstly, aerodynamic theory is used to create formulae which can be used to calculate both the axial and rotational energy losses in the wake of a propeller. These are then validated with wind tunnel data, and it is shown that, for an optimal propeller, the loss in efficiency due to rotational velocity is always very small. However, for a propeller with high-blade settings, this efficiency loss can become quite high. The load distribution along a propeller blade is also examined, and the maximum thrust generation is found to always be between 0.7R-0.85R (where R is the radius of the blade). Finally, the benefits of counter-rotating propellers over single propellers with an increased diameter are examined, and it was found that only small gains in propeller efficiency can be expected with counter-rotating blades.

A book by Theodorsen (1948) entitled Theory of Propellers [12] was a development of the earlier propeller theory developed by Drzewiecki (1892), and Reynolds et al. (1890). Theodorsen focussed on the wake downstream of the propeller more than the propeller itself, as “thrust, torque and efficiency are uniquely and completely established by the knowledge of the wake far behind the propeller and of the wake alone. To establish these quantities, in fact, it is absolutely unnecessary to know the design of the propeller.” Theodorsen introduces a mass coefficient, which represents the effective cross-section of the column of air behind the propeller divided by the projected propeller wake area. Formulae are then developed for calculating the torque, thrust, energy loss and efficiency, in terms of the mass coefficient, for any propeller with ideal circulation distribution, from the downstream wake conditions.

In 1979, McMasters and Henderson working for Boeing, published a report entitled Low-Speed Single-Element Airfoil Synthesis [13] which outlines the design criteria and the practical restraints for an airfoil operating at low Reynolds numbers. The report also describes a ‘synthesis’ based approach to airfoil design, in contrast to the typical ‘analytical’ approach. First, the boundary layer characteristics are defined, as they define the pressure distribution and therefore the airfoil performance. From this, powerful computational tools were used to generate an airfoil shape which yields the pre-defined boundary layer.

A Master’s thesis by Gamble (2009) entitled Automatic Dynamic Propeller Testing at Low Reynolds Numbers [14] focuses on the testing of propellers and reducing propeller test times, by designing a system in which propeller data can be acquired systematically. The effect of the Reynolds number on efficiency and propeller pitch were also examined.

Another Master’s thesis by Koch (1998) called the Design and Performance Calculations of a Propeller for Very High-Altitude Flight [15] examined propellers for use in low Reynolds number flow regimes (60,000-100,000). A design process utilizing BET and ADPAC (a numerical Navier-Stokes code) was used to optimise a propeller for high altitude flight. The operating conditions for this study were 25.9 km (85,000 ft) at Mach 0.4, and each propeller was required to absorb 63.4 kW (85 bhp).

Propeller Performance Measurement for Low Reynolds Number Unmanned Aerial Vehicle Applications [16] focused on the testing of small propellers, particularly those used on UAV systems. An Integrated Propulsion Test System (IPTS) was used to evaluate the performance of many propellers, and a database of performance data was created. This IPTS setup is similar to the one in use at the University of Illinois at Urbana-Champaign’s (UIUC) Department of Aerospace Engineering, where a vast propeller database has been created [2]. Using a load cell to measure thrust and a torque cell to measure torque, efficiency curves have been plotted, and in modern propeller research this database is often used for comparative reasons.

Similarly, in Low Reynolds Propellers for Increased Quadcopters Endurance [17], the design process for creating a propeller for use in low Reynolds number flow is detailed by Carvalho (2013). A validation exercise of the UIUC Propeller database [2] was also completed, by comparing the UIUC experimental results of the APC Slow Flyer 10” x 7” to simulated results on JBLADE software. The highest efficiency value was found to be 61% from JBLADE simulations, compared to the experimental value of 59%. Note: The UIUC database was updated in March 2015, to account for a new wind tunnel correction method for the drag due to the propeller mounting fixture. This had the result of increasing the peak efficiency of props by approximately 10%.

A journal paper entitled Propeller Performance Measurements at Low Reynolds Numbers [18] by Silvestre, et al. (2015) discussed the design of propellers for use on high-altitude airships. Two APC propellers were tested between 4 m/s and 28 m/s in 1 m/s intervals, and the data was compared to the UIUC database in order to validate the test-rig. It was shown that as the Reynolds number increases with increasing RPM values, the efficiency and thrust coefficient of the propeller was increased significantly. It is also stated that the test-rig would be used to create experimental data that could improve the JBLADE software in the future.

In a 2016 paper by Kuantama, et al. entitled Quadcopter Propeller Design and Performance Analysis [19] used SolidWorks to analyse the efficiency of a two-bladed 16” x 5” propeller. The range of RPM values used in this case study was 1,000-9,000 RPM. The experimental data was compared to datasheets of existing propellers in order to validate the SolidWorks flow simulation method.

A book by Prior (2018) Optimising Small Multi-Rotor Unmanned Aircraft—A Practical Design Guide [20] gives an overview of propellers used on small multi-rotor aircraft. Propeller variables are discussed, including the geometry and airfoil choice that needs to be made when designing a propeller. Basic propeller theory, including momentum theory, various propeller coefficients, and propeller efficiency is outlined. The effect of using overlapping rotors is explored, from large aircraft such as the Boeing CH-47 Chinook helicopter, to smaller aircraft such as the Malloy Aeronautics Hoverbike. Finally, Prior discusses the challenges of propeller design in low Reynolds number flow, using the analogy that at low Reynolds numbers, air becomes more viscous and therefore achieving high efficiencies is much harder.

Finally, an Experimental Study on the Aerodynamic and Aeroacoustic Performances of a Bio-Inspired UAV Propeller [21] discussed the performance of a baseline propeller and a bio-inspired propeller with the same design thrust (3N) and solidity (0.12) are analysed. It was found that for the same thrust generation, the bio-inspired propeller emitted lower noise under constant power input. The planform shape of the bio-inspired propeller was inspired by the cicada wing and the shape of maple seed.

1.6. The Optimum Propeller

There is ample literature describing and validating various methods of analysing propeller performance, however there is significantly less existing literature on design methods for optimum propellers.

All propellers have spanwise or radial circulation distributions which minimise the kinetic energy loss associated with the generation of lift (thrust). By determining the geometry which results in these circulation distributions, one can determine the geometry of the propeller which results in the minimum losses (most efficient). In 1979, Larrabee published Design of Propellers for Motorsoarers [22] which explains his method, and forms a series of equations which are solved iteratively to form the final propeller geometry. The two geometric parameters this method calculates are the chord distribution and the twist distribution. The suitability of the method for three very different applications was then discussed. These were propellers for use on a human-powered aircraft, a powered hang glider and a motorsoarer. Comparisons were made between the calculated propellers and the actual propellers used on these aircraft, and in all instances the calculated geometry was found to be very close to the actual optimum geometry.

The Design of Optimum Propellers (Adkins & Liebeck, 1994) [23] further explains the process, which is based around the Betz condition that a constant-displacement velocity across the propeller wake provides a propeller design with maximum efficiency. Using the circulation functions developed by Betz & Prandtl [24] and Goldstein [25], one can use the Betz condition to obtain the chord and twist distributions of a theoretical propeller with minimum losses. The original equations stated by Betz and Glauert [10] were only applicable for propellers with light-loading, however, Theodorsen [12] showed that the Betz condition for minimum loss can be applied to heavy-loading as well, by including the contraction of the wake in the equations.

In his paper entitled Increasing of Aircraft Propeller Efficiency by Using Variable Twist Propeller Blades (Klesa, 2008) [26], the method is simplified to a simple step by step iterative process. Klesa shows that a MATLAB code could be written to generate the optimum chord and twist distributions, and shows that the varying flight conditions will present different optimal geometry, meaning that a propeller can only be optimised for one design point at a time, unless deploying variable pitch mechanisms.

2. Materials and Methods

2.1. Propeller Theory—Geometry

2.1.1. Airfoil Geometry

An airfoil is a 2D cross-sectional shape of an object that, when moved through a fluid, produces an aerodynamic force. The most common application of airfoils is on the wings of fixed-wing aircraft, however they are also used on propellers. Among others, the main objective when selecting an airfoil is to create the highest lift to drag ratio. The most suitable airfoil will depend on the flow conditions the airfoil will be operating in. Using online tools such as JavaFoil [27] or QProp [28], existing airfoils can be analysed, and new airfoils can be created and analysed. Figure 3 depicts a NACA0012 in JavaFoil, which is one of the most common airfoils. As can be seen, the airfoil is symmetric (it has no camber). This can be deduced from the ’00’ digits in the airfoil name. The airfoil has a thickness to chord ratio of 12% (hence the ’12’ digits in the airfoil name).

However, most airfoils used on propellers are cambered, as symmetric airfoils require either high angles of attack or high flow speeds to produce significant lift, resulting in poor lift to drag ratios. The NACA3414 airfoil, as seen in Figure 4, is a cambered airfoil. It has 3% maximum camber (signified by the ’3’ digit), with the maximum camber occurring at 40% chord (signified by the ’4’ digit), and a thickness to chord ratio of 14%.

2.1.2. Propeller Geometry

Aside from the diameter, the pitch of the propeller is the most important geometric parameter when designing a propeller. The pitch of a propeller is defined as the distance that, if the propeller were moving through a soft solid, the propeller would move forward in one complete rotation. Of course, in conventional propellers the pitch will decrease from the hub to the tip, and so it has become convention to specify the pitch of a fixed pitch propeller at 75% of the distance from the hub to the blade tip. Propellers with high pitch values can produce higher thrust levels at higher speeds than low pitch propellers, whilst low pitch propellers can produce higher thrust levels at low speeds. Therefore, the main driver in selecting the pitch for a propeller is the velocity of the propeller. From the pitch, the theoretical angle of twist can be calculated for the length of the blade using Eqn. 2.

$$Angle\mathrm{ }of\mathrm{ }Twist=ta{n}^{-1}\left(\frac{p}{2\mathsf{\pi }r}\right)$$

(2)

Another important geometric aspect of a propeller is the taper profile. Whilst rotor blades found on helicopters tend to have no taper, and wings on fixed wing aircraft tend to have linear tapers, the taper profile of propellers is typically quite complicated, involving complex curves. The taper profile is defined as the variation of chord length along the blade. Typically, the chord length is small at the hub, it increases, resulting in a maximum chord length between 25% and 75% along the blade radius, and then it decreases towards the blade tip.

2.2. Non-Dimensional Coefficients

The performance parameters of a propeller (or rotor) can be expressed in non-dimensional form. Understanding how to calculate these quantities is therefore key to understanding propeller theory. It should be noted that these coefficients differ for propellers and rotors. A rotor is defined as an aerodynamic propulsion device which operates in static conditions, whereas a propeller operates in non-static conditions. Simply put, for the calculation of these coefficients a rotor is a propeller operating with a forward flow speed of zero. From Leishman (2008) [29], the non-dimensional coefficients can be found using Eqns. 3 to 8, and are categorised in Table 1 below.

$$C_T=\frac{T}{\frac{1}{2}\mathsf{\rho }A{\left(\mathsf{\Omega }R\right)}^2}$$

(3)

$$C_{Tp}=3.875C_T=\frac{T}{\mathsf{\rho }{n}^2{D}^4}$$

(4)

$$C_P=\frac{P}{\frac{1}{2}\mathsf{\rho }A{\left(\mathsf{\Omega }R\right)}^3}$$

(5)

$$C_{Pp}=12.175C_P=\frac{P}{\mathsf{\rho }{n}^3{D}^5}$$

(6)

$$C_Q=\frac{Q}{\frac{1}{2}\mathsf{\rho }A{\left(\mathsf{\Omega }R\right)}^{2}R}$$

(7)

$$C_{Qp}=1.938 C_Q=\frac{Q}{\mathsf{\rho } n^2{D}^5}$$

(8)

Another dimensionless quantity of rotor blades is the blade solidity, which can be defined as the ratio of total blade area to the disk area, as defined in Eqn. 9. Values of the solidity for conventional helicopters vary from about 0.06 to 0.12.

$$\mathsf{\sigma }=\frac{Area of rotor blades}{Area of rotor disk}=\frac{BcR}{\mathsf{\pi } R^2}=\frac{Bc}{\mathsf{\pi } R}$$

(9)

2.3. Propeller/Rotor Efficiency

The efficiency calculation of a blade propulsion system varies depending on whether the system is a rotor or propeller driven system. Rotor driven systems, such as those found on helicopters, use a hover efficiency metric called the Figure of Merit (FoM), whilst propeller driven systems use a more conventional propulsive efficiency term η_P.

Figure of Merit

The Figure of Merit is a dimensionless factor used to calculate hover efficiency of rotor driven propulsion systems. The concept of a Figure of Merit was first proposed by Renard (1903) and then Glauert (1935), but was first presented in its present form by Richard H. Prewitt of Kellett Aircraft. The Figure of Merit is defined in Eqn. 10.

$$FoM=\frac{Ideal Power Required to hover}{Actual Power Required to hover}=\frac{C_T^{\frac{3}{2}}}{\sqrt{2}{C}_P}$$

(10)

Eqn. 10 is derived from simple momentum theory, and thus assumes there is no viscous losses when calculating the ideal power required to hover. In the actual power required to hover measurement, viscous effects will result in both profile and induced loss contributions, and therefore for a real rotor, the Figure of Merit must be less than 1. Hover propulsion systems can be optimised to achieve FoM values of up to 0.9, however, achieving an optimal FoM value will result in a significantly lower propulsive efficiency value. Therefore, rotor propulsion systems are often designed to have FoM values of around 0.6, resulting in higher efficiencies at high speeds at the expense of lower efficiencies in hover conditions.

The disc loading of a rotor is defined as the ratio of the thrust produced by the rotor to the rotor disk area. To make a meaningful comparison between the FoM values of multiple rotors, the rotors must have the same disc loading. The disc loading will increase the induced power relative to the profile power, resulting in a higher FoM which is a misleading value for comparison. Generally, rotors with high disc loadings will give higher figure of merit values than rotors with low disc loadings.

Propulsive Efficiency

The propulsive efficiency, η_P, of a propeller driven propulsion system is a more useful expression of propeller performance, since in their operation, all propulsion systems will have axial flow with a non-zero velocity. Firstly, the advance ratio, J which is needed to calculate the propulsive efficiency, is defined in Eqn. 11.

$$J=\frac{Axial Flow Velocity}{Tip Rotational Velocity}=\frac{V}{nD}$$

(11)

Typically, in the range of 0-2, the advance ratio is the ratio of axial velocity to rotational velocity at the tip of the propeller. As can be seen in Eqn. 12, the propulsive efficiency is directly proportional to the advance ratio, J.

$${\mathsf{\eta }}_P=J\frac{C_T}{C_P} \equiv \frac{Tc}{\mathrm{P}\mathrm{c}}$$

(12)

As can be seen in Eqn. 12, the propulsive efficiency of a propeller is directly affected by the coefficients of thrust and power. Looking back to Eqns. 4 and 6, it can be seen that both the angular velocity and the propeller diameter will affect the propulsive efficiency. The fluid density, which is in the denominator of both coefficients, cancels out and therefore does not affect the efficiency. However, the propulsive efficiency is determined by other factors not explicitly seen in Eqns. 4 and 6. The propeller pitch/diameter ratio, airfoil geometry, spanwise chord, taper profile and material choice all affect the thrust and power terms seen in Eqns. 4 and 6; thus, they all affect the propulsive efficiency. The efficiency curve of an APC 10” x 7” propeller is shown in Figure 5, with data from the UIUC propeller database [2].

Well-designed propellers typically have maximum efficiencies at advance ratios of between 0.5 and 0.7. However, the maximum efficiency for a propeller only occurs at one advance ratio (J), meaning that as the advance ratio increases, the efficiency will quickly reduce to a minimum. This sudden drop means that designing an efficient propeller that will be operated over a range of advance ratios is very challenging, as one can select an advance ratio at which the propeller will operate at maximum efficiency, but any variation in advance ratio can lead to a drastic drop in efficiency. Therefore, the maximum efficiency of a propeller is limited to a single airspeed at a given angular velocity. Propellers on larger aircraft suffer from compressibility losses, whereby the tips of the blades approach the speed of sound, resulting in shock wave formation and the losses associated with these. Transonic flow begins at about Mach 0.8.

2.4. Reynolds Number

The Reynolds number of a certain flow condition is defined as the ratio of inertial forces to viscous forces, and it is a dimensionless parameter that can help predict flow patterns. The Reynolds number is defined in Eqn. 13, where L is the characteristic length.

$$Re=\frac{\mathsf{\rho }\cdot V\cdot L}{\mathsf{\mu }}$$

(13)

As can be seen in Eqn. 13, if all other parameters are held constant, the Reynolds number will decrease if the characteristic length decreases. When analysing propeller flow, the characteristic length is the chord length of the blades, and so smaller blades will result in smaller Reynolds numbers. As can be seen in Figure 6, at a Reynolds number of 100,000 for smooth airfoils, there is a sudden transition from a low maximum section C_L/C_D ratio to a high maximum section C_L/C_D ratio. Propellers used on small drones often operate in Reynolds numbers less than 100,000, and the maximum lift/drag ratio that can be achieved is therefore limited by their small diameters.

When operating below the 100,000 Reynolds number threshold, the lower lift/drag ratios associated with these flow regimes, are the biggest contributor to the lower efficiencies associated with scale propellers. The velocity component in Eqn. 13 means that the angular velocity of the propeller affects the Reynolds number. The Reynolds number will vary along the length of the blade, due to varying chord lengths and rotational velocities, and at low angular velocities sometimes no part of the blade will exceed the 100,000 Reynolds number threshold. It is therefore critical for a propellers efficiency to operate at the highest RPM value possible (assuming tip speeds do not approach Mach 0.8).

2.5. Propeller Geometry Optimization

The following section describes a method that is used to iteratively design propeller chord and twist distributions. The following systems of equations are only reviewed briefly in this paper, but a more detailed explanation and derivation has been written by both Larrabee (1979) [22] and Adkins & Liebeck (1994) [23]. The method below assumes that the shaft power supplied to the motor is known, however, if the required thrust of the propeller is known instead, Klesa (2008) [26] describes a method that is largely the same as the one described here. This method assumes light-loading (<400 N/m²) and incompressible flow.

2.6. Calculating Propeller Geometry

The numerical method is based on the Betz criteria [24]. The Betz criterion for optimum loading for a lightly loaded actuator disk in axial flight states that the trailing edge vortex sheets, after an initial deformation close to the rotor disk, move axially downward as a rigid screw surface (i.e., the ratio of translational velocity to rotational velocity is constant). For every iteration of this method, the displacement velocity is assumed to be radially constant. If the wake helix angle is φ_s, the axial and swirl velocity components are therefore given by Eqns. 14 and 15.

$$w_{axial}=v^{\prime }{\mathit{Cos}}^{2}2{\mathsf{\varphi }}_s$$

(14)

$$w_{swirl}=v^{\prime }\cos {\mathsf{\varphi }}_s\sin {\mathsf{\varphi }}_s$$

(15)

In this method, the displacement velocity will be represented by the displacement velocity ratio, which is defined in Eqn. 16. This is the variable which will be optimised to result in optimum propeller geometry, but this will be discussed later.

$$\mathsf{\zeta }=\frac{v^{\prime }}{V}$$

(16)

Both Betz & Prandtl (1919) and Goldstein (1929) derived momentum loss functions, which describe the energy lost in a propeller wake. Whilst Goldstein’s function is more accurate, Prandtl’s momentum loss function, F, will be used here for simplicity.

$$F\left(r\right)=\frac{2}{\mathsf{\pi }}\arccos e^{-f}$$

(17)

$$f\left(r\right)=\frac{B}{2}\cdot \frac{1-\mathsf{\xi }}{\sin {\mathsf{\varphi }}_t}$$

(18)

$$\mathsf{\xi }\left(r\right)=\frac{r}{R}$$

(19)

ξ is the non-dimensional radius and is calculated by dividing the radial co-ordinate, r, by the propeller tip radius, R as in Eqn. 19. φ_t is the flow angle at the propeller tip, and is a function of the displacement velocity, as seen in Eqn. 20. From the flow tip angle, the flow angle, φ, of every radial point can be found using Eqn. 21.

$$\tan {\mathsf{\varphi }}_t=\frac{V}{\mathsf{\Omega }R}\cdot \left(1+\frac{\mathsf{\zeta }}{2}\right)$$

(20)

$$\mathsf{\varphi }\left(r\right)=\arctan \frac{\tan {\mathsf{\varphi }}_t}{\mathsf{\xi }}$$

(21)

From Eqns. 17-21, Eqns. 22 and 23 can be derived, which describe the circulation, Γ, at any radial point r. The circulation function, G, is more accurate for vortex sheets which are flat, parallel and closely spaced, however, it is an adequate approximation for most propellers.

$$\mathsf{\Gamma }\left(r\right)=\frac{2\mathsf{\pi }{V}^2\mathsf{\zeta }G}{B\mathsf{\Omega }}$$

(22)

$$G\left(r\right)=F\cdot \frac{\mathsf{\Omega }r}{V}\cdot \cos \mathsf{\varphi }\sin \mathsf{\varphi }$$

(23)

Applying the Kutta-Joukowski theorem for a blade element at radius r, the lift per unit length of the blade can be equated as in Eqn. 24, where W is the local total velocity.

$$\frac{\mathsf{\rho }{W}^{2}c{c}_L}{2}=\mathsf{\rho }W$$

(24)

Applying the expression for Γ in Eqn. 22, results in Eqn. 25.

$$Wc\left(r\right)=\frac{4\mathsf{\pi }{V}^{2}G\mathsf{\zeta }}{c_{L}B\mathsf{\Omega }}$$

(25)

The total velocity W can be expressed by Eqn. 26.

$$W\left(r\right)=\frac{V\left(1+a\right)}{\sin \mathsf{\varphi }}$$

(26)

a and a^, are variables which appear several times in this method. The derivation of them will not be included in this paper, but can be found in ’Increasing of Aircraft Propeller Efficiency by Using Variable Twist Propeller Blades’, a paper by Klesa (2008) [26].

$$a\left(r\right)=\frac{\mathsf{\zeta }}{2}\cdot {cos}^2\mathsf{\varphi }\cdot \left(1-\mathsf{ϵ}\tan \mathsf{\varphi }\right)$$

(27)

$$a^{\prime }\left(r\right)=\frac{\mathsf{\zeta }}{2}\cdot \frac{V}{\mathsf{\Omega }r}\cdot \cos \mathsf{\varphi }\sin \mathsf{\varphi }\left(1+\frac{\mathsf{ϵ}}{\tan \mathsf{\varphi }}\right)$$

(28)

E is the section drag to lift ratio and is the only airfoil property that is required for this method. By dividing Eqn. 25 by Eqn. 26, the value of c can be calculated for any radial co-ordinate. Using Eqn. 29, the angle of twist (β) can be found at any radial co-ordinate.

$$\mathsf{\beta }\left(r\right)=\mathsf{\varphi }+\alpha$$

(29)

2.7. Optimising Propeller Geometry

The two main geometric parameters of the propeller, the chord and twist distribution, have now been calculated, however these parameters are not optimised. The following section details how the displacement velocity ratio is optimised to produce optimal propeller geometry in an iterative process. Firstly, the non-dimensional design power coefficient, P_c, has to be calculated using Eqn. 30, where P is the shaft power supplied to the propeller.

$$P_c=\frac{2P}{\mathsf{\rho }{V}^3\mathsf{\pi }{R}^2}$$

(30)

The key equation in the iterative process is Eqn. 31, which after simple integration results in Eqn. 32.

$$\frac{d{P}_c}{d\xi }=\frac{d{J}_1}{d\xi }\cdot \mathsf{\zeta }+\frac{d{J}_2}{d\xi }\cdot {\mathsf{\zeta }}^2$$

(31)

$$P_c=J_1\cdot \mathsf{\zeta }+J_2\cdot {\mathsf{\zeta }}^2$$

(32)

J₁ and J₂ are simply variables used in the iteration process and have no physical meaning. ’Design of Propellers for Motorsoarers’, a paper by Larrabee (1979) [22], contains a full derivation of Eqns. 31, 32, 33 and 34.

$$\frac{d{J}_1}{d\mathsf{\xi }}=4\mathsf{\xi }G\left(1+\frac{\mathsf{ϵ}}{\tan \mathsf{\varphi }}\right)$$

(33)

$$\frac{d{J}_2}{d\mathsf{\xi }}=\frac{1}{2}\cdot \frac{d{J}_1}{d\mathsf{\xi }}\cdot \left(1-\mathsf{ϵ}\tan \mathsf{\varphi }\right)\cdot {\mathrm{c}\mathrm{o}\mathrm{s}}^2\mathsf{\varphi }$$

(34)

After evaluating the derivatives of J₁ and J₂ at every radial co-ordinate, integrating with respect to the non-dimensional radius will result in values for J₁ and J₂. Re-arranging Eqn. 32, results in Eqn. 35. Using calculated values of P_c, J₁ and J₂, a new value of the displacement velocity ratio can be found.

$$\zeta =-\frac{J_1}{2J_2}+\sqrt{{\left(\frac{J_1}{2J_2}\right)}^2+\frac{P_c}{J_2}}$$

(35)

The process is then repeated using the new value for the displacement velocity ratio, and is iterated until the displacement velocity ratio converges. Initially, an estimate for the displacement velocity ratio is used, which has to be small and is often zero. If the thrust, T is known then a slightly different equation is used with derivatives of I₁ and I_2. in the above procedure [22,23].

3. Results

3.1. Implementing the Optimization Method Using MATLAB

The required inputs for a MATLAB computer program which implements the iterative method are found in Table 2, along with their units. The non-dimensional hub radius is found by dividing the hub radius by the propeller tip radius.

The angle of attack input specifies the additional twist angle that is desired on the propeller. As can be seen in Eqn. 29, if the angle of attack is specified as zero, the radial twist angle distribution will be the same as the radial flow angle distribution.

As the method is not analytical, there will inevitably always be errors. The main sources of error in the method are the assumptions of light-loading and incompressible flow, made throughout the method. Without greatly increasing the method complexity, there is little that can be done about this. However, the accuracy of the method can be increased by varying two of the inputs in Table 2. Firstly, the number of blade elements defines the blade element width. Using more blade elements will result in smaller blade element widths, resulting in a more complete and accurate output, as can be seen in Figure 7a and Figure 7b. The displacement velocity ratio convergence level defines how close the calculated value of zeta needs to be to the original value for the program to classify the displacement velocity ratio as ’converged’.

Using a smaller value for this input will therefore increase accuracy, at the inevitable expense of increased computation time, due to an increased number of iterations, as can be seen in Figure 8.

3.2. Reynolds Number Distribution

Using Eqn. 13, the Reynolds number distribution along the blade can be calculated, using the optimal chord distribution output from the program. The re-arrangement of Eqn. 13 to an equation suitable for the computer program is found below.

$$Re\left(r\right)=\frac{\mathsf{\rho }\cdot W\left(r\right)\cdot c\left(r\right)}{\mathsf{\mu }}=\frac{\mathsf{\rho }\cdot \mathsf{\Omega }r\cdot c\left(r\right)}{\mathsf{\mu }}$$

(36)

As already discussed in this paper, a Reynolds number of 100,000 represents a transition point at which the lift-to-drag ratio of a smooth airfoil increases drastically. In propeller design, it is therefore highly desirable for the Reynolds number to be above 100,000 at most radial co-ordinates. Whilst this could be achieved at all radial co-ordinates by simply increasing the chord length close to the propeller hub, structural requirements dictate that this chord profile is not commonly seen in small (less than 21-inch diameter) propellers. There would also be a drag penalty associated with increasing the chord length, which would reduce the lift-to-drag ratio of the blade sections close to the hub.

’Increasing of Aircraft Propeller Efficiency by Using Variable Twist Propeller Blades’, a paper by Klesa (2008) [26], details a method which, given that the above radial variables have been calculated at every radial co-ordinate, can be easily used to calculate the radial distribution of thrust generated by the propeller. After derivation, the equation used in the computer program is shown in Eqn. 37.

$$dT\left(r\right)=2\mathsf{\pi }\mathsf{\rho }V\cdot \left(1+a\right)\cdot \left(2VaF\right)\cdot dr$$

(37)

An example thrust distribution, plotted by calculating dT at every radial co-ordinate, is shown in Figure 9.

As can be seen in Figure 9, the majority of the thrust developed by a propeller is between the non-dimensional radii of 0.4 and 0.95. Using numerical integration, it is found that for the example in Figure 9, 85% of the thrust is developed between 0.4 and 0.95 of the propeller radii, but other propellers will yield very similar results. Therefore, in the computer program created for this paper, the program ensures that all blade stations between non-dimensional radii of 0.4 and 0.95 experience Reynolds numbers above the 100,000 threshold.

The unmodified Reynolds number distributions of two example propellers are found in Figure 10.

For clarification, a flow chart outlining the basic architecture of the MATLAB program which outputs optimal propeller geometry and subsequently modifies this to minimise the effect of low Reynolds number on propeller efficiency is shown in Figure 11.

3.3. Program User Interface Design

In order to make the MATLAB program user friendly, MATLAB libraries inputdlg() and listdlg() were used to create a series of pop-up windows that ask the user for the inputs required. The user can choose between whether the required thrust or consumed shaft power is known.

Following this, the pop-up window in Figure 12 appears, asking the user for the required inputs. Lift and drag coefficient data can be found in ‘Aerodynamic data banks for Clark-Y, NACA 4-digit and NACA 16-digit airfoil families’ [29].

Following this, a window appears asking the user if they require the Reynolds number chord adjustment. If the user does require this, they are then prompted to input the desired minimum Reynolds number, which by default is 100,000. If not, the code is run as normal. The MATLAB code can be found via a link in Appendix A.

3.4. Validation of the MATLAB Program

In order to validate the MATLAB program, its output was compared to online data for the APC ’Thin Electric’ 10” x 7”, from the UIUC propeller database [2]. Inputs to the program, found in Table 3, were chosen to match the performance of the APC ’Thin Electric’ 10” x 7” propeller. For example, the propeller RPM and forward velocity were chosen to match the advance ratio at which the APC propeller has a maximum efficiency (advance ratio of 0.6). The airfoil coefficients of lift and drag are from the Clark-Y airfoil, which is a common airfoil for propeller design [30,31]. An angle of attack of zero was used, as due to assumptions made in the method the twist angles are large, thus eliminating the requirement for additional angle of attack. The shaft power was approximated by comparing the required shaft powers of propellers of a similar diameter and pitch on the APC database [32].

In Figure 13, Figure 14 and Figure 15, the chord, twist and Reynolds number distribution for both the MATLAB program output and the APC ’Thin Electric’ 10” x 7” propeller are plotted.

In Table 4, calculated parameters for both the MATLAB output and the UIUC database wind tunnel measurements can be found. The C_T and C_P values in Table 4 are for an advance ratio of 0.5751, which is the advance ratio at which the maximum efficiency occurs. Results from the APC propeller database [32] are also included for reference, however these are simulated and not experimental.

As can be seen in Figure 13, Figure 14 and Figure 15, the MATLAB program output closely resembles the geometry of the APC propeller. The APC propeller has larger chord lengths towards the propeller tip, and this is most likely due to structural issues, which the MATLAB program does not consider. Finally, there is a significant difference (16%) between the maximum propulsive efficiencies, the main source of which being the assumptions of light-loading, incompressible flow and zero tip losses made in the MATLAB program method. On the contrary, there is only a 2% difference between the UIUC wind tunnel test data and the APC simulation database. The smaller difference is most likely due to similar assumptions made in both calculation processes. However, the fact that the MATLAB program computes an optimized geometric profile in terms of the blade twist and blade chord, as well as computes these using a baseline Re number of 100,000 may mean that this is a better propeller design in terms of propulsion efficiency. Note that the MATLAB code produces a higher pitch over the spanwise radius, r of the propeller.

4. Discussion

In this paper, the efficiency of propellers used for small drones has been examined, and various ways of improving their efficiency have been evaluated.

This was first done by developing a MATLAB code, based on an optimisation method described in various literature on propeller theory. For a set of inputs specified by the user, the code will output a full chord and twist distribution for an optimal propeller blade. A function allowing the user to modify the chord distribution to ensure the blade was operating above a given Reynolds number threshold was then added. Another function was programmed which enabled the user to view the thrust and torque distributions along the blade. Finally, a case study comparing results from the MATLAB code to published results for a similar propeller, the APC ’Thin Electric’ 10” x 7”. This case study validated the MATLAB code, proving that it could be a useful tool to use for designing optimum propellers.

Summary:

• A MATLAB tool has been developed, which takes inputs from the end user and calculates the optimal propeller geometry. The code is both fast, reliable and easy to use.
• Theoretically, the propeller design created will achieve higher efficiencies at a range of RPM values than commercially available propellers. From the MATLAB code, it was calculated that whilst operating at 6,500 RPM, an efficiency of 0.87 can be expected.

4.1. Recommendations for Further Work

The research and design work that has been undertaken for this paper has highlighted several areas in which further research could be beneficial.

The MATLAB tool developed in this study outputs a plot of chord and twist distributions for an optimum design based on end user input. These distributions then have to be manually converted into a set of airfoil curves, which are then linked by a surface loft in the SolidWorks CAD program. If a program which automated this process, and automatically output a SolidWorks file could be created, this would reduce the design time of optimal propellers to a fraction of the time it currently takes to complete the manual process. The MATLAB code could also be improved by incorporating tip loss calculations into the optimisation process, as in the current code tip losses are neglected. Furthermore, the code could also be designed to calculate the optimal blade winglet profile to reduce the losses from wingtip vortices. Finally, the code could be modified so that the user is warned when the propeller tip speed is approaching transonic speeds of Mach 0.8. This would make it easier for the user to avoid the compressibility effects and shock waves that are associated with such tip speeds.