Wind Energy Systems
BEM - Exercise 3
This is a small notebook that contains the calculation and visualisation for one of the exercise sheets in the course Wind Energy Systems 2021. All of the derivation was done by hand. This is notebook only serves the purpose of validating the derivations and sharing the results.
# Handling all of the imports here
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import math
from scipy.integrate import quad
mpl.style.use('ggplot')
print('Matplotlib version: ', mpl.__version__)
Matplotlib version: 3.4.2
#######################################
# Defining functions to reduce clutter
#######################################
# Defining the functional axial induction factor as it also changes with mu
def a_axial(mu):
return ((0.8 + (28 * mu)) * (10**(-2)))
# Defining the functional tangential induction factor as it also changes with mu
def a_tang(mu):
return ((0.3 + (0.6 / mu) + (2.8 * mu)) * (10**(-3)))
# Defining the functional solidity as it also obviously changes with mu
def solidity(mu):
return (8 / 441 * mu)
# Calculating the effective wind magnitutde based on mu
def W(mu, tip_speed_ratio):
return (12 ** 2) * (((mu**2) * (tip_speed_ratio**2) * (1 + a_tang(mu))**2) + (1-a_axial(mu))**2)
########################################
# Defining constants
########################################
B = 3 # Number of Blades in Turbine A
tip_speed_ratio = 7 # The tip speed ratio for Turbine A
R = 50 # The radius of the rotor blade
phi = 5 # Flow Angle or Effective Velocity Angle
c_l = 1 # The 2D lift coefficient
c_d = 0.01 # The 2D drag coefficient
u_inf = 12 # The free stream wind velocity
row = 1.225 # The density of free stream air
########################################
# Defining the evaluation function
########################################
# Setting up mu
# Have to remove the 0 to avoid division by it in tangential induction factor calculation
step_size = 100
mu = np.linspace(0, 1, step_size)
mu = mu[1:len(mu)]
dmu = 1 / step_size
# Calculating wind speed vs mu just for visualisation
w = W(mu, tip_speed_ratio)
# Showing Effective Wind Speed vs mu
plt.figure(figsize = (21, 3))
plt.title('\nEffective Wind Speed Magnitude (squared) vs. $\mu$\n')
plt.xlabel('$\mu$')
plt.ylabel('W - Effective Wind Speed Magnitutde ($\\frac{m}{s}$)')
plt.plot(mu, w)
plt.show()
# Okay now that we have setup the functions and constants. Now comes the main event
def thrust_coefficient(mu):
return (4 * a_axial(mu) * (1 - a_axial(mu)))
def thrust(mu, B, u_inf, row, phi, c_l, c_d, R, tip_speed_ratio):
# First we calculate the constants
constants = (R ** 2) * (3.14) * row * (u_inf ** 2) * (1/B) * ((c_l * math.cos(phi)) + (c_d * math.sin(phi)))
# Then we combine everything we know
return (constants * W(mu, tip_speed_ratio) * solidity(mu) * mu)
dT = thrust(mu, B, u_inf, row, phi, c_l, c_d, R, tip_speed_ratio) * dmu
plt.figure(figsize = (21, 3))
plt.title("\nThrust vs. $\mu$\n")
plt.xlabel('$\mu$')
plt.ylabel('T - Thrust Magnitude (N)')
plt.plot(mu, dT)
plt.show()
########################################
# Thrust and Thrust Coefficient
########################################
# Evaluating the integrals
thrust_one_blade, err = quad(thrust, 0, 1, args=(B, u_inf, row, phi, c_l, c_d, R, tip_speed_ratio))
C_T, err = quad(thrust_coefficient, 0, 1)
print('Calculation Results:')
print()
print('Thrust on whole rotor of turbine A:\t\t\t' + str(thrust_one_blade * 3) + ' N')
print('Thrust Coefficient on whole rotor of turbine A:\t\t' + str(C_T))
print()
Calculation Results:
Thrust on whole rotor of turbine A: 9984620.280029813 N
Thrust Coefficient on whole rotor of turbine A: 0.47825066666666666