# Rotor Blade Element Analysis (RBEA)

A project log for Inexpensive Composite Propellers/Rotors

How to make inexpensive propellers using mostly hardware store materials and a CNC router

Peter McCloud 02/24/2017 at 04:200 Comments

Recently, the rotor designs were re-visited to double check some of the design assumptions. In the process the code was polished up and some improvements were made.

Bug Fixes

Found that the code would not converge for some cases and the user wouldn't know. Not so much a bug as poor programming. Added logic to notify the user if the case doesn't converge

Improvements

Broke out the inputs into a file (inputs.txt), keeping the code cleaner. Additionally switched to using splines for the airfoil data instead of fits done in a spreadsheet. Added support for Mach and alpha airfoil coeffs vs just alpha.

## Python Code (repo location)

```import sys
import math
import numpy as np
import scipy.interpolate
# Python Script for determining Rotor Performance using Blade Element Theory
# Definitions
# theta - The physical angle the airfoil is rotated from the plane of rotation
# phi	- The change in flow angle due to the induced airflow.  Makes the angle of attack smaller
# alpha_rad	- The effective angle of attack the airfoil experiences
# alpha_rad + phi = theta
pi = math.pi
execfile('inputs.txt')
# Determine airfoil method
try:
airfoil_name
except:
# airfoil_name is not defined, see if tip & root airfoil names exist
try:
tip_airfoil_name
except NameError:
print 'Neither airfoil_name nor tip_airfoil_name in inputs file'
print 'Exiting!'
sys.exit(1)
else:
# tip_airfoil_name exists, see if root_airfoil_name exists
try:
root_airfoil_name
except:
print 'root_airfoil_name is not defined in the inputs file'
print 'Exiting!'
sys.exit(1)
else:
airfoil_type = 'blended'
else:
airfoil_type = 'single'
if airfoil_type == 'single':
# Create splines
mach_data = airfoil[:,0]
alpha_data = airfoil[:,1]
cd_data = airfoil[:,2]
cl_data = airfoil[:,3]
cm_data = airfoil[:,4]
# Determine Mach Range
if min(mach_data) == max(mach_data):
# Single Mach Number present, Create 1D splines
airfoil_data_type = '1D'
cd_spline = scipy.interpolate.splrep(alpha_data,cd_data)
cl_spline = scipy.interpolate.splrep(alpha_data,cl_data)
cm_spline = scipy.interpolate.splrep(alpha_data,cm_data)
else:
# Multiple Mach Numbers present, Create 2D splines
airfoil_data_type = '2D'
cd_spline = scipy.interpolate.bisplrep(mach_data,alpha_data,cd_data)
cl_spline = scipy.interpolate.bisplrep(mach_data,alpha_data,cl_data)
cm_spline = scipy.interpolate.bisplrep(mach_data,alpha_data,cm_data)
else:
# Tip
# Create splines
mach_data = airfoil[:,0]
alpha_data = airfoil[:,1]
cd_data = airfoil[:,2]
cl_data = airfoil[:,3]
cm_data = airfoil[:,4]
# Determine Mach Range
if min(mach_data) == max(mach_data):
# Single Mach Number present, Create 1D splines
tip_airfoil_data_type = '1D'
tip_cd_spline = scipy.interpolate.splrep(alpha_data,cd_data)
tip_cl_spline = scipy.interpolate.splrep(alpha_data,cl_data)
tip_cm_spline = scipy.interpolate.splrep(alpha_data,cm_data)
else:
# Multiple Mach Numbers present, Create 2D splines
tip_airfoil_data_type = '2D'
tip_cd_spline = scipy.interpolate.bisplrep(mach_data,alpha_data,cd_data)
tip_cl_spline = scipy.interpolate.bisplrep(mach_data,alpha_data,cl_data)
tip_cm_spline = scipy.interpolate.bisplrep(mach_data,alpha_data,cm_data)
# Root
# Create splines
mach_data = airfoil[:,0]
alpha_data = airfoil[:,1]
cd_data = airfoil[:,2]
cl_data = airfoil[:,3]
cm_data = airfoil[:,4]
# Determine Mach Range
if min(mach_data) == max(mach_data):
# Single Mach Number present, Create 1D splines
root_airfoil_data_type = '1D'
root_cd_spline = scipy.interpolate.splrep(alpha_data,cd_data)
root_cl_spline = scipy.interpolate.splrep(alpha_data,cl_data)
root_cm_spline = scipy.interpolate.splrep(alpha_data,cm_data)
else:
# Multiple Mach Numbers present, Create 2D splines
root_airfoil_data_type = '2D'
root_cd_spline = scipy.interpolate.bisplrep(mach_data,alpha_data,cd_data)
root_cl_spline = scipy.interpolate.bisplrep(mach_data,alpha_data,cl_data)
root_cm_spline = scipy.interpolate.bisplrep(mach_data,alpha_data,cm_data)
# Initial Calculations
area = pi*diameter*diameter/4
avg_chord = (root_chord + tip_chord)/2
theta_inc = (tip_theta-root_theta)/(num_elements-1)
chord_inc = (tip_chord-root_chord)/(num_elements-1)
# Initialize lists
theta = []
chord = []
thrust = []
drag =[]
torque = []
# Initialize total_output file
total_output = open('./output.txt','w')
total_output.write(' RPM\tT(N)\tP(kW)\tM(N-m)\\(kg-cm)\t\tT(lbf)\tP(Hp)\tM(ft-lb)\t\tFM\tM_tip\tv_induced\n')
# Initialize total_output file
for i in range(num_elements):
# Populate radii and theta lists
for i in range(num_elements):
chord.append(root_chord+chord_inc)
theta.append(root_theta+i*theta_inc)
# Specified Theta
#theta = [18.00, 17.34, 16.68, 16.03, 15.37, 14.71, 14.05, 13.39, 12.74, 12.08, 11.42, 10.76, 10.10, 9.44, 8.79, 8.13, 7.47, 6.82, 6.16, 5.50]
#for i in range(num_elements):
#print theta
for RPM in RPMs:
# Initial Calcs at RPM
n = RPM/60					# RPS
omega = n*2*pi				# angular velocity
V_tip = omega*tip_radius	# tip velocity
# Initialize/Clear Lists
thrust = []
drag = []
moment = []
torque = []
for i in range(num_elements):
# Guess at initial values of inflow and swirl factor
# Note: swirl currently isn't used
v_inflow = 10
#v_swirl = 0.1
# Iterate at Each Blade Element
iter = 0
finished = False
while( finished == False):
v_axial = v_inflow				# axial velocity
mach = v_mag/a
alpha = min(alpha,12)
# Find section coefficients
if airfoil_type == 'single':
if airfoil_data_type == '1D':
cd = scipy.interpolate.splev(alpha,cd_spline)
cl = scipy.interpolate.splev(alpha,cl_spline)
cm = scipy.interpolate.splev(alpha,cm_spline)
elif airfoil_data_type == '2D':
cd = scipy.interpolate.bisplev(mach,alpha,cd_spline)
cl = scipy.interpolate.bisplev(mach,alpha,cl_spline)
cm = scipy.interpolate.bisplev(mach,alpha,cm_spline)
else:
# Tip
if tip_airfoil_data_type == '1D':
tip_cd = scipy.interpolate.splev(alpha,tip_cd_spline)
tip_cl = scipy.interpolate.splev(alpha,tip_cl_spline)
tip_cm = scipy.interpolate.splev(alpha,tip_cm_spline)
elif tip_airfoil_data_type == '2D':
tip_cd = scipy.interpolate.bisplev(mach,alpha,tip_cd_spline)
tip_cl = scipy.interpolate.bisplev(mach,alpha,tip_cl_spline)
tip_cm = scipy.interpolate.bisplev(mach,alpha,tip_cm_spline)
# Root
if root_airfoil_data_type == '1D':
root_cd = scipy.interpolate.splev(alpha,root_cd_spline)
root_cl = scipy.interpolate.splev(alpha,root_cl_spline)
root_cm = scipy.interpolate.splev(alpha,root_cm_spline)
elif root_airfoil_data_type == '2D':
root_cd = scipy.interpolate.bisplev(mach,alpha,root_cd_spline)
root_cl = scipy.interpolate.bisplev(mach,alpha,root_cl_spline)
root_cm = scipy.interpolate.bisplev(mach,alpha,root_cm_spline)
# Interpolate cd, cl, cm
tip_portion = float(i)/(num_elements-1)
root_portion = float(num_elements-1-i)/(num_elements-1)
cd = tip_portion*tip_cd + root_portion*root_cd
cl = tip_portion*tip_cl + root_portion*root_cl
cm = tip_portion*tip_cm + root_portion*root_cm
# Appply drag factor
cd = cd*drag_factor
q = 0.5*rho*pow(v_mag,2)								# local dynamic pressure
# momentum check on inflow and swirl factors
# increment iteration count
iter += 1
# check for convergence
if ( math.fabs(v_inflow_new-v_inflow)< 0.001):
finished = True
# check to see if iterations are stuck
elif(iter>maximum_iterations):
finished=True
print 'RPM:',RPM,'element:',i,'exceed maximum number of iterations (',maximum_iterations,')'
v_inflow = v_inflow + relaxation_factor*(v_inflow_new-v_inflow)
#v_swirl = v_inflow + 0.5*(v_swirl_new-v_swirl)
thrust.append(DtDr*r_inc)
drag.append(DdDr*r_inc)
moment.append(DmDr*r_inc)
torque.append(DqDr*r_inc)
#print RPM,i,v_inflow
for i in range(num_elements):
theta[i] = theta[i]*180/pi
# Totals
T = sum(thrust)			# total thrust (N)
P = sum(torque)*omega/1000	# total power (kW)
M = sum(moment)			# total moment (N-m)
M_kg = M*10.1971621298		# total moment (kg-cm)
# Ideals
v_ideal = pow((T/(2*rho*pi*pow(tip_radius,2))),0.5) # Ideal Induced Velocity (m/s)
P_ideal = T*v_ideal/1000 # Ideal Power (kW)
# Compute Coefficients
ct = T/(rho*n*n*pow(diameter,4))
cq = sum(torque)/(rho*n*n*pow(diameter,5))
# Adjust for Tip Loss using Empirical Equations
B = 1 - pow(2*ct,2)/2
P = P_ideal/B + (P - P_ideal)
T_imp = T*0.224808943			# N to lbf
P_imp = P*0.00134102209*1000	# kW to Hp
M_imp = M*0.737562149			# N-m to ft-lbf
FM = P_ideal/P
M_tip = v_mag/a
total_output.write('{:5d}\t{:5.2f}\t{:5.2f}\t{:6.3f}\\{:6.3f}\t\t{:5.2f}\t{:5.2f}\t{:6.3f}\t\t{:5.3f}\t{:5.3f}\t{:5.3f}\n'.format(RPM,T,P,M,M_kg,T_imp,P_imp,M_imp,FM,M_tip,v_ideal))
# Close total_output
total_output.close()