Code Updates

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.

Additions

Code for calculating blade moments

Python Code (repo location)

import sys
import math
import numpy as np
import scipy.interpolate
# Rotor Blade Element Analysis
# 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
# Read inputs file
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' 
# Read airfoil data
if airfoil_type == 'single':
 airfoil = np.loadtxt('./airfoils/'+airfoil_name+'.dat',skiprows=1,delimiter=',')
 # 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
 airfoil = np.loadtxt('./airfoils/'+tip_airfoil_name+'.dat',skiprows=1,delimiter=',')
 # 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 
 airfoil = np.loadtxt('./airfoils/'+root_airfoil_name+'.dat',skiprows=1,delimiter=',')
 # 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
tip_radius = diameter/2
area = pi*diameter*diameter/4
avg_chord = (root_chord + tip_chord)/2
sigma=2*avg_chord/2/pi/tip_radius			# solidity
# Blade Element Calculation
r_inc = (tip_radius-root_radius)/(num_elements-1)
theta_inc = (tip_theta-root_theta)/(num_elements-1)
chord_inc = (tip_chord-root_chord)/(num_elements-1)
# Initialize lists
radius = []
theta = []
theta_rad = []
alpha_rad = []
chord = []
thrust = []
drag =[]
torque = []
# Initialize total_output file
total_output = open('./output.txt','w')
# Print Data Header
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
radial_output = open('./radial_output.txt','w')
# Print Data Header
radial_output.write('i\\RPM\t')
for i in range(num_elements):
 radial_output.write(str(i)+'\t')
radial_output.write('\n')
# Populate radii and theta lists
for i in range(num_elements):
 radius.append(root_radius+i*r_inc)
 chord.append(root_chord+chord_inc)
 theta.append(root_theta+i*theta_inc)
 theta_rad.append(theta[i]/180*pi)
 alpha_rad.append(0)
# 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):
# theta_rad[i] = (theta[i]/180*pi)
#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
   v_radial = omega*radius[i]			# disk plane velocity
   #v_radial = omega*radius[i] - v_swirl	# disk plane velocity
   v_mag = pow((v_axial*v_axial+v_radial*v_radial),0.5)	# local velocity at blade (Resultant of Induced & Radial)
   mach = v_mag/a
   phi_rad = math.atan2(v_axial,v_radial)	# flow angle (radians)
   alpha_rad[i] = theta_rad[i]- phi_rad		# blade angle of attack
   #cl,cd = airfoil_coef(alpha_rad[i])
   alpha = alpha_rad[i]*180/pi
   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
   DtDr = q*blades*chord[i]*(cl*math.cos(phi_rad)-cd*math.sin(phi_rad))			# thrust grading
   DdDr = q*blades*chord[i]*(cd*math.cos(phi_rad)+cl*math.sin(phi_rad))			# drag grading
   DmDr = q*blades*chord[i]*cm													# moment grading
   DqDr = q*blades*chord[i]*radius[i]*(cd*math.cos(phi_rad)+cl*math.sin(phi_rad))	# torque grading
   # momentum check on inflow and swirl factors
   v_inflow_new = DtDr/(4*pi*radius[i]*rho*v_axial)
   #v_swirl_new = DqDr/(4*pi*pow(radius[i],3)*rho*v_axial*omega)
   # 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,')'
   # Updates Values
   v_inflow = v_inflow + relaxation_factor*(v_inflow_new-v_inflow)
   #v_swirl = v_inflow + 0.5*(v_swirl_new-v_swirl)    
  phi = phi_rad*180/pi
  alpha = alpha_rad[i]*180/pi 
  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
 # Convert radians to degrees
 radial_output.write('{:d}\t'.format(RPM))
 for i in range(num_elements):
  theta[i] = theta[i]*180/pi
  alpha_rad[i] = alpha_rad[i]*180/pi
  radial_output.write('{:5.3f}\t'.format(alpha_rad[i]))
 radial_output.write('\n')
 # 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()
radial_output.close()