Preliminary sizing of a rotor blade#
As previously mentioned we will use the work of Adkins and Liebeck to perform both optimisation and off design analysis. Let us start with the rotor blade optimisation.
The method has been implemented in Python in the package called tuduam, i.e TU Delft Urban Air Mobility. Hence, when creating the input we have to keep in mind the documentation on what the functions expect. Feel free to roam around the documentation to see what other functionalities are offered by the package. Keep in mind however that the package is under active development.
The design is initiated with the specified conditions of power (or thrust), hub and tip radius, rotational rate, freestream velocity, number of blades, and a finite number of stations at which blade geometry is to be determined. Let us define this in the next block
Defining the conditions#
Most of the code below does not need explaining, however as previously mentioned we must take into account the API documentation . Here it says that the blade parameter must be defined in the Propeller data structure. Furthermore, BEM analysis requires airfoil polar at a large range of reynolds numbers. tuduam relies on the user to input this data in the form of a path to a directory which contains all the polars. These files should be in the common format outputted by xfoil. A tutorial on how to create this data has also been made which can be found in the course. Sample data for the NACA 4412 is readily available however in the tuduam repository. The NACA 4412 is shown in Figure 1
Figure 1: The NACA 4412 airfoil
The conditions required for the optimisation are loaded into the Python kernel below. Carefully read the comments surrounding the code for more pointers. Note that the repository makes use of data structures in order to limit input size into functions, in this case we have to make use of the Propeller data structure.
!pip install tuduam
!git clone https://github.com/saullocastro/tuduam.git
# If the following error occurs "fatal: destination path 'tuduam' already exists and is not an empty directory." please continue as your environment is already set up correctly
Requirement already satisfied: tuduam in c:\users\damie\anaconda3\envs\tuduam_notebook\lib\site-packages (2024.15)
Requirement already satisfied: numpy in c:\users\damie\anaconda3\envs\tuduam_notebook\lib\site-packages (from tuduam) (1.26.4)
Requirement already satisfied: scipy in c:\users\damie\anaconda3\envs\tuduam_notebook\lib\site-packages (from tuduam) (1.12.0)
Cloning into 'tuduam'...
import os
import tuduam as tud
import rich
import tuduam.propulsion as prop
# Define parameters required according to API documentation
data_path = os.path.realpath(os.path.join(os.path.abspath('.'), 'tuduam', 'examples', 'data')) # Path to airfoil data
coord_path = os.path.realpath(os.path.join(os.path.abspath('.'), 'tuduam', 'examples', 'naca_4412.txt')) # Path to airfoil coordinates
n_stations = 25 # number of stations
# Define operational parameters
thrust_req = 400 # N
v_cruise = 210/3.6 #m/s
h_cruise = 1200 # m
isa = tud.ISA(h_cruise) # ISA altitude class
a_snd = isa.soundspeed() # Speed of sound m/s
rho_cr = isa.density() # Air density kg/m^3
dyn_vis = isa.viscosity_dyn() # Dynamic viscosity Pa*s
# Define the specified conditions for the propeller geometry in the data structure
prop_class = tud.Propeller(
n_blades=6, # Number of blades
r_prop=0.55, # propeller radius m
rpm_cruise=1200, # revolutions per minute
tc_ratio=0.12, # Thickness to chord ratio of airfoil
xi_0=0.1 # nondimensional hub radius (r_hub/R)
)
rich.print(f"{a_snd=} \n{rho_cr=} \n{dyn_vis=} \n{prop_class=}")
a_snd=335.62572904948746 rho_cr=1.0899402822825992 dyn_vis=1.7736932645281385e-05 prop_class=Propeller(n_blades=6.0, r_prop=0.55, rpm_cruise=1200.0, xi_0=0.1, chord_arr=None, rad_arr=None, pitch_arr=None, tc_ratio=0.12)
Running the optimisation#
Having defined the conditions we now run the actual optimisation. In the code it will be as simple as calling a single line however there is much going on in the background. A general overview of procedure is described here however for a full description refer to the paper (download link available in the course in sources).
The general procedure can be condensed to the following:
Select an initial estimate for \(\zeta\), the displacement velocity ratio (\(\zeta\) = 0 will work).
Determine the values for \(F\), the Prandtl momentum loss factor, and \(\phi\), the flow angle, at each blade station.
Determine the product \(Wc\), the local velocity multiplied by the local chord. Additionally, compute the Reynolds number.
Determine \(\epsilon\), the drag-to-lift ratio, and \(\alpha\), the angle of attack, from airfoil section data.
If \(\epsilon\) is to be minimized, change \(C_l\), and repeat Steps 3 and 4 until this is accomplished at each station.
Determine \(\alpha\) and \(\alpha'\), the axial and rotational interference factor, respectively.
Compute the chord from Step 3, and the blade twist \(\beta = \alpha + \phi\).
Determine the four derivatives in I and J, and numerically integrate these from the non-dimensional hub radius, \(\xi_0\) to the non-dimensional radius \(\xi\) = 1.
Determine \(\zeta\) and \(T_c\), the thrust coefficient.
If this new value for \(\zeta\) is not sufficiently close to the old one (e.g., within 0.1%), start over at step 2 using the new \(\zeta\).
Determine propeller efficiency as \(\eta = \frac{T}{P_c}\), and other features such as solidity.
Let us now run the optimisation and plot the resulting blade geometry. The plotting is performed using a class called PlotBlade, which has a method for both a 3D and 2D plot. Feel free to play with the parameters of the rotor by rerunning the above code box with new parameters. Note that the optimisation will crash however for infeasible designs.
import numpy as np
# Initialise BEM class
design = prop.BEM(data_path, prop_class, rho_cr,
dyn_vis, v_cruise, n_stations, a_snd, T=thrust_req)
# Execute the BEM optimisation with initial estimate of 0
res_dict = design.optimise_blade(zeta_init=0)
# Print the resulting dictionary
rich.print(res_dict)
rich.print(f'power required = {np.round(res_dict["pc"]*rho_cr*(prop_class.rpm_cruise/60)**3*32*prop_class.r_prop**5/1e3,2)} Kw')
# Plot the resulting blade geometry
plot = prop.PlotBlade(prop_class, coord_path)
{ 'chord_arr': array([0.06450542, 0.0815267 , 0.0970076 , 0.11083296, 0.12292547, 0.13324371, 0.14178268, 0.1484006 , 0.15344058, 0.15682568, 0.15863253, 0.1589442 , 0.15784599, 0.15542109, 0.15174425, 0.14687961, 0.14087312, 0.13374916, 0.12549776, 0.11606204, 0.1054556 , 0.0931046 , 0.07862803, 0.06075668, 0.03493964]), 'pitch_arr': array([1.50066203, 1.46814551, 1.43597542, 1.4042232 , 1.37295481, 1.34048482, 1.31210253, 1.27738256, 1.24858189, 1.22049744, 1.19315549, 1.16657602, 1.14077307, 1.11575523, 1.09152613, 1.06808489, 1.0436814 , 1.02179809, 1.00067843, 0.98030877, 0.96416404, 0.94524556, 0.92702529, 0.90948357, 0.89260002]), 'alpha_arr': array([0.04363323, 0.04537856, 0.04712389, 0.04886922, 0.05061455, 0.05061455, 0.05410521, 0.05061455, 0.05235988, 0.05410521, 0.05585054, 0.05759587, 0.05934119, 0.06108652, 0.06283185, 0.06457718, 0.06457718, 0.06632251, 0.06806784, 0.06981317, 0.07504916, 0.07679449, 0.07853982, 0.08028515, 0.08203047]), 'station_arr': array([0.0649, 0.0847, 0.1045, 0.1243, 0.1441, 0.1639, 0.1837, 0.2035, 0.2233, 0.2431, 0.2629, 0.2827, 0.3025, 0.3223, 0.3421, 0.3619, 0.3817, 0.4015, 0.4213, 0.4411, 0.4609, 0.4807, 0.5005, 0.5203, 0.5401]), 'drag_to_lift_arr': array([0.02467022, 0.02453238, 0.02441314, 0.02412773, 0.02379947, 0.02349407, 0.02362806, 0.00833238, 0.00831749, 0.0082426 , 0.00819419, 0.00815943, 0.00812609, 0.00811742, 0.00807504, 0.00803422, 0.00794982, 0.00791248, 0.00788762, 0.00788581, 0.02165054, 0.02161504, 0.02146528, 0.02128869, 0.02137814]), 'v_e': 84.42019742325951, 'solidity': array([0.94912391, 0.91915302, 0.88646351, 0.85146964, 0.81460912, 0.77631708, 0.73703042, 0.69637412, 0.65617987, 0.61603246, 0.57619974, 0.53689612, 0.49828701, 0.46049087, 0.42357523, 0.38756479, 0.35243364, 0.31810968, 0.28445652, 0.25126068, 0.21849139, 0.18495599, 0.15001845, 0.11150943, 0.06177541]), 'cl': array([0.74153846, 0.75207692, 0.76261538, 0.77315385, 0.78369231, 0.79423077, 0.80476923, 0.81530769, 0.82584615, 0.83638462, 0.84692308, 0.85746154, 0.868 , 0.87853846, 0.88907692, 0.89961538, 0.91015385, 0.92069231, 0.93123077, 0.94176923, 0.95230769, 0.96284615, 0.97338462, 0.98392308, 0.99446154]), 'cd': array([0.01829392, 0.01845023, 0.01861784, 0.01865444, 0.01865146, 0.01865971, 0.01901514, 0.00679345, 0.00686897, 0.00689399, 0.00693985, 0.0069964 , 0.00705345, 0.00713147, 0.00717933, 0.0072277 , 0.00723556, 0.00728496, 0.00734519, 0.00742662, 0.02061798, 0.02081196, 0.02089398, 0.02094644, 0.02125974]), 'eff': 0.7873405743678134, 'tc': 0.2269752189628377, 'pc': 0.28828086136051734, 'zeta': 0.44719577704754987 }
power required = 4.05 Kw
#Plot 2D outputs
plot.plot_blade() # Creates two 2D plots
#Plot 3D outputs
plot.plot_3D_plotly() # Creates an interactive 3D plot
Offdesign analysis#
We now have an optimised blade which is stored in the Propeller data structure and in the output of the function, i.e res_dict. We can do a quick test in order to check this. Run the following cell to see what the data structure now contains. As seen it now contains 3 extra arrays describing the propeller blade.
We would now like to test whether it can meet a certain off design condition. If you want to know more about a specific step, please refer to the paper which is available on the course site. The general procedure for doing so is as follows:
Specify a freestream velocity, RPM and \(\Delta \beta\)
Make in initial estimate on the flow angle using \(\phi = \frac{\lambda}{\xi}\) where \(\xi\) is the nondimensional spanwise location and \(\lambda= \frac{V_{inf}}{\omega*R} \)
Since \(\beta\) is known, compute \(\alpha = \beta - \phi\) and the Reynolds number at each station.
The airfoil coefficients are then known from the airfoil section data.
The new estimate for \(\phi\) can then be computed as in equation 1. Let a and \(a`\) be the axial and rotational interference factor, r the local radius and V the freestream velocity. Note that a and \(a`\) are clipped to -0.7 and 0.7 respectively. This is required for adequate convergence near the tip (Viterna and Janetzke (1981)). For more details please see Adkins and Liebeck (1984) on the course site.
Steps 1 to 5 are then repeated until convergence is reached.
The thrust and torque per unit radius can then be found using:
Then finally the thrust and power coefficient are defined as folows. Where n is the revolutions per second.
Finally, this can then be integrated. In the source code this is done using a simple trapezoid integration, more specifically the scipy.integrate.trapezoid function.
Hover example#
Let’s do an sample example of a vtol rising vertically to see how this would look like using the tuduam repository. The inflow from the freestream is now significantly reduced. The parameter that significantly changes is the freestream velocity, the only significant component left is the induced velocity caused by the rotor. Determining the value of the induced velocity is difficult, engineering judgement is used. [1] mentions an equation to estimate it however it is only reliable for higher flight speeds. Thus in our script we have assumed the freestream velocity to be 10 m/s.
Furthermore it is important to consider how we should change our pitch control of the rotor. As seen in Figure 2, the resultant relative wind will lay nearly in parrallel to tip-path axis. Therefore if we would keep the same pitch as the one during cruise the blade will most surely stall. The implication in the code being an error that an angle of attack outside of the dataset has been encountered.
Therefore we would like to lower our pitch angle. To increase the thrust for hover we would like to increase our RPM. However keep in mind that our rpm is limited by the onset of shock waves at the tip of the blade. Thus we have to compute the maximum rpm and limit the rpm to that value. We will now analyse the propeller at a singular offdesign conditions. Namely a freestream velocity of 10 m/s, a maximum mach number of 0.7 at the tip and a 39 degree decrease in pitch on the rotor blades. Read and run the code below for the results.
Sources#
[1] Castles, W. Jr., & De Leeuw, J. H. (1954). The normal component of the induced velocity in the vicinity of a lifting rotor and some examples of its application (Report No. 1184). Georgia Institute of Technology.
import numpy as np
from scipy.optimize import minimize
# Define constants required for analysis
mach_max = 0.7 # Point when shockwave occurs over the airfoil
isa = tud.ISA(0) # Assume sea level
rho = isa.density()
dyn_vis = isa.viscosity_dyn()
a_snd = isa.soundspeed()
# Operations conditions
vtip_max = mach_max*a_snd
omega_max = vtip_max/prop_class.r_prop
rpm_hvr = omega_max/(2*np.pi)*60 # Vertical speed of the vtol
v_hvr = 10 # Vertical speed of the vtol plus additional
delta_pitch = np.radians(-39)
rich.print(f"{rpm_hvr=}")
hover_cond = prop.OffDesignAnalysisBEM(data_path, prop_class, v_hvr, rpm_hvr, rho, dyn_vis, a_snd)
res = hover_cond.analyse_propeller(delta_pitch, abs_extrapolation=7)
power_req = res["power_coeff"]*rho*(rpm_hvr/60)**3*32*prop_class.r_prop**5/1e3
rich.print(res)
rich.print(f"The power required for a single rotor is {power_req} kW")
rpm_hvr=4135.43320654644
{ 'thrust': 4880.664989838316, 'torque': 294.0218860337935, 'eff': 0.3833096128070276, 'thrust_coeff': 0.5728390932248073, 'power_coeff': 0.1971153712449012, 'AoA': array([0.12277787, 0.21520819, 0.27347206, 0.31892348, 0.34060297, 0.34950426, 0.34547205, 0.35195982, 0.33585566, 0.32746425, 0.31953733, 0.30883281, 0.29703362, 0.28429357, 0.27132936, 0.25596598, 0.24239548, 0.22354675, 0.2082226 , 0.19300554, 0.18230526, 0.17562534, 0.16223733, 0.14912745, 0.13591189]), 'lift_coeff': array([1.19482265, 1.35056012, 1.40478863, 1.1920388 , 1.19503635, 1.19855278, 1.55793305, 1.55460521, 1.58987287, 1.6089598 , 1.63197297, 1.6621659 , 1.70106215, 1.72257239, 1.73685942, 1.75631737, 1.76430709, 1.73880235, 1.71841164, 1.6939743 , 1.67804967, 1.74210487, 1.73295633, 1.7226558 , 1.67964973]), 'drag_coeff': array([0.02334077, 0.04634767, 0.08521674, 0.14794629, 0.14831832, 0.14875475, 0.11732507, 0.02192104, 0.10682134, 0.02643343, 0.09170785, 0.08135113, 0.06912503, 0.06115158, 0.05001949, 0.0637648 , 0.081345 , 0.10155172, 0.02521023, 0.022538 , 0.18440827, 0.02069385, 0.01885703, 0.01682833, 0.03141791]) }
The power required for a single rotor is 127.32957449453333 kW
Finding the max thrust of a blade#
Let us now try a range of values for the both the change in pitch and rpm values to see if we can achieve better thrust values. We could opt to use scipy.optimize for this purpose however the analysis is prone to failure due to either to high AoA or a non-converging solution. This will complicate the optimisation where we’ll have to assign penalization when such an error occurs. Instead we’ll use a brute force method since the analysis is swift. We’ll create a mesh where x will be the rpm and y the pitch. Then for each coordinate we’ll compute the thrust, after which a surface plot can be created. This is implemented below, the cell takes about 30 seconds to run, so some patience is required.
Feel free to play with the parameters!
import plotly.graph_objects as go
n = 14 # Mesh size
pitch_arr = np.linspace(-20, -50, n)
rpm_arr = np.linspace(2500,rpm_hvr, n)
X,Y = np.meshgrid(rpm_arr,pitch_arr)
Z_thrust = np.zeros(X.shape)
Z_eff = np.zeros(X.shape)
Z_power = np.zeros(X.shape)
for col, rpm in enumerate(rpm_arr):
hover_cond = prop.OffDesignAnalysisBEM(data_path, prop_class, v_hvr, rpm, rho, dyn_vis, a_snd)
for row, dpitch in enumerate(np.radians(pitch_arr)):
try:
res = hover_cond.analyse_propeller(dpitch, abs_extrapolation=7)
except ValueError:
continue
except RuntimeError as error:
continue
power_req = res["power_coeff"]*rho*(rpm_hvr/60)**3*32*prop_class.r_prop**5/1e3
Z_thrust[row,col] = res["thrust"]
Z_eff[row,col] = res["eff"]*100
Z_power[row,col] = power_req
# Create the figure
fig = go.Figure(data=[go.Surface(z=Z_thrust, x=X, y=Y)])
fig.update_layout(
scene=dict(
xaxis_title='RPM',
yaxis_title= 'Delta pitch [deg]',
zaxis_title=f'Thrust [N]'))
fig.update_layout(
scene_camera=dict(
eye=dict(x=1.8, y=1.8, z=1.8)))
fig.update_layout(title=f'Thrust plot {n} x {n} mesh ', autosize=True,
margin=dict(l=65, r=50, b=65, t=90))
# Show the plot
fig.show()
Iteration of the rotor design.#
So what happens if the maximum thrust we find is not sufficient? Remember when a thrust or power value was specified for which the rotor would be optimized? We can increase this value to increase the maximum thrust value, this however does come at a cost that the blade is less efficient at the cruise setting. With the updated geometry we can then perform the offdesign analysis again to find the newly acquired maximum thrust. This process is then repeated until the engineer is sufficiently happy with the margin on the maximum thrust and the design thrust of the propeller.
The opportunity is given below to go through this process yourself with the code cell below where all of the relevant code has been appended. The mesh spacing has been increased to decrease computation time, however when nearing to your final design you can decrease the mesh spacing again. The design point is described in the first few lines of code. Sometimes the results are odd with the majority failing or having no solution at all. In this case it is advised to increase the radius of your propeller or change the amount of blades to reduce disk and/or blade loading.
# Define operatonal parameters
thrust_req = 350
v_cruise = 210/3.6
h_cruise = 1200
isa = tud.ISA(h_cruise)
a_snd = isa.soundspeed() # Speed of sound
rho_cr = isa.density() # Air density
dyn_vis = isa.viscosity_dyn() # Dynamic viscosity
# Define the specified conditions for the propeller geometry in the data structure
prop_class = tud.Propeller(
n_blades=6, # Number of blades
r_prop=0.525, # propeller radius
rpm_cruise=1200,
tc_ratio=0.12, # Thickness to chord ratio of airfoil
xi_0=0.1 # nondimensional hub radius (r_hub/R)
)
#---------------------------- Blade geometry optimization ---------------------------------
# Initiazlie BEM class
design = prop.BEM(data_path, prop_class, rho_cr,
dyn_vis, v_cruise, n_stations, a_snd, T=thrust_req)
# Execute the BEM optimisation with initial estimate of 0
res_dict = design.optimise_blade(zeta_init=0)
plot = prop.PlotBlade(prop_class, coord_path)
plot.plot_3D_plotly() # Creates an interactive 3D plot
#------------------- Off design analysis ---------------------------------------------------------
n = 5 # Mesh size
pitch_arr = np.linspace(-20, -50, n)
rpm_arr = np.linspace(2500,rpm_hvr, n)
X,Y = np.meshgrid(rpm_arr,pitch_arr)
Z_thrust = np.zeros(X.shape)
Z_eff = np.zeros(X.shape)
Z_power = np.zeros(X.shape)
for col, rpm in enumerate(rpm_arr):
hover_cond = prop.OffDesignAnalysisBEM(data_path, prop_class, v_hvr, rpm, rho, dyn_vis, a_snd)
for row, dpitch in enumerate(np.radians(pitch_arr)):
try:
res = hover_cond.analyse_propeller(dpitch, abs_extrapolation=7)
except ValueError:
continue
except RuntimeError as error:
continue
power_req = res["power_coeff"]*rho*(rpm_hvr/60)**3*32*prop_class.r_prop**5/1e3
Z_thrust[row,col] = res["thrust"]
Z_eff[row,col] = res["eff"]*100
Z_power[row,col] = power_req
# Create the figure
fig = go.Figure(data=[go.Surface(z=Z_thrust, x=X, y=Y)])
fig.update_layout(
scene=dict(
xaxis_title='RPM',
yaxis_title= 'Delta pitch [deg]',
zaxis_title=f'Thrust [N]'))
fig.update_layout(
scene_camera=dict(
eye=dict(x=1.8, y=1.8, z=1.8)))
fig.update_layout(title=f'Thrust plot {n} x {n} mesh ', autosize=True,
margin=dict(l=65, r=50, b=65, t=90))
# Show the plot
fig.show()