import pandas as pd import math import numpy as np import matplotlib.pyplot as plt from scipy import stats import statistics import matplotlib.pyplot as plt from scipy.interpolate import * from scipy import integrate as intg

df = pd.read_csv("/work/DU25_A17.dat",error_bad_lines=False, sep=" ",header=None) df.columns = ['alpha', 'Cl', 'Cd', 'moment', 'NaN'] df.head()

#donnée pour le pitch (beta) et la vitesse de rotation en fonction de la vitesse du vent v_rot_rads = ((2*np.pi)/60)*np.array([7,7.2,7.5,7.9,8.5,9.2,10.35, 11.5, 11.9, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1, 12.1]) data = {'Vitesse du vent' : [3,4,5,6,7,8,9, 10, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0], 'pitch' : [0,0,0,0,0,0,0,0, 0.00, 3.83, 6.60, 8.70, 10.45, 12.06, 13.54, 14.92, 16.23, 17.47, 18.70 ,19.94 ,21.18 ,22.35 ,23.47], 'vitesse de rotation' : v_rot_rads } df_pitch = pd.DataFrame(data)

# Données de longueur de corde en fonction du rayon data_bis = {'Rayon' : [2.8667, 5.6000, 8.3333, 11.7500, 15.8500, 19.9500, 24.0500, 28.1500, 32.2500, 36.3500, 40.4500, 44.5500, 48.6500, 52.7500, 56.1667, 58.9000, 61.6333], 'Corde' : [2.7333, 2.7333, 2.7333, 4.1000, 4.1000, 4.1000, 4.1000, 4.1000, 4.1000, 4.1000, 4.1000, 4.1000, 4.1000, 4.1000, 2.7333, 2.7333, 2.7333], 'Twist' : [13.308, 13.308, 13.308, 13.308, 11.480, 10.162, 9.011, 7.795, 6.544, 5.361, 4.188, 3.125, 2.319, 1.526, 0.863, 0.370, 0.106]} df_corde = pd.DataFrame(data_bis)

#Calcul de a et a' vent = np.array(['trois','quatre','cinq','six','sept','huit','neuf','dix','onze','douze','treize','quatorze','quinze','seize','dixsept','dixhuit','dixneuf','vingt','vingtun','vingtdeux','vingttrois','vingtquatre','vingtcinq']) B = 3 rho = 1.225 for i in range (len(data['Vitesse du vent'])): U = data['Vitesse du vent'][i] b=data['pitch'][i]+data_bis['Twist'][i] w = data['vitesse de rotation'][i] exec('{} = pd.DataFrame()'.format(vent[i])) #permet de former une nouvelle dataframe pour chaque vitesse du vent avec comme nom ceux qui se trouvent dans le vecteur vent globals()[vent[i]] = (globals()[vent[i]]).assign(r=[0],a=[0],a_bis=[0],CL=[0],CD=[0],Tforce = [0], Qforce = [0], Qmoment = [0], Tmoment=[0], SCE=[0])#le globals permet de faire en sorte d'appeler les différents dataframes produits lors des différentes boucles sur base d'un string se trouvant dans un vecteur. globals()[vent[i]]=(globals()[vent[i]]).drop((globals()[vent[i]]).iloc[0,:])#on supprime la première ligne de 0 qui ne sert à rien (ça a juste permis de placer les colonne dans la dataframe) for j in range (len(data_bis['Rayon'])): r = data_bis['Rayon'][-j-1] c = data_bis['Corde'][-j-1] dfa = pd.DataFrame(columns=['a','a_bis','CL','CD','Tforce','Qforce','Qmoment','Tmoment','SCE']) if j == 0: #lors de la première itération on ne peut pas se baser sur les valeurs précédentes for q in range (20,101,1): dfa_bis = pd.DataFrame(columns=['a_bis','CL','CD','Tforce','Qforce','Qmoment','Tmoment','SCE'])#je le fais en deux fois parce que a et a' ne vont pas toujours varier sur la même gamme de valeur ni sur le même nombre de valeurs for q_bis in range (1,q+1,1): #a_bis n'est pas censé être plus grand que a try: a = q/100 #obligé de faire ça comme ça: la boucle fonctionne pas si je veux la faire tourner avec des nombres à virgule a_bis = q_bis/100 phi = math.degrees(math.atan((U*(1-a))/(w*r*(1+a_bis)))) cosphi = math.cos(math.radians(phi)) sinphi = math.sin(math.radians(phi)) #v_res = ((1-a)*U)/(sinphi) v_res = (((1-a)*U)**2)/(w*r*(1+a_bis))**2 alpha = phi - b smallest_gap = 3000 for x in df['alpha']: #boucle pour trouver la valeur du tableau la plus proche du alpha trouvé if abs(alpha - x)<smallest_gap: smallest_gap = abs(alpha - x) alpha_tableau = x CL = df[df['alpha'] == alpha_tableau]['Cl'].iloc[0]#on prend la valeur de CL pour l'alpha correspondant. Le [0] permet d'avoir un chiffre et pas une dataframe. CD = df[df['alpha'] == alpha_tableau]['Cd'].iloc[0] Tforce = 0.5*B*rho*(v_res**2)*(CL*cosphi+CD*sinphi)*c*r Qforce = 0.5*B*rho*(v_res**2)*(CL*sinphi-CD*cosphi)*c*((r**2)/2) Tmoment = 4*np.pi*rho*(U**2)*(1-a)*a*((r**2)/2) Qmoment = 4*np.pi*rho*U*(1-a)*a_bis*w*((r**4)/4) SCE = (((Tforce-Tmoment)**2)+((Qforce-Qmoment)**2))**(0.5) dfa_bis_new = pd.DataFrame(data=np.array([[a_bis,CL,CD,Tforce,Qforce,Qmoment,Tmoment,SCE]]), columns=['a_bis','CL','CD','Tforce','Qforce','Qmoment','Tmoment','SCE']) dfa_bis = pd.concat([dfa_bis,dfa_bis_new], ignore_index=True) #on ajoute une nouvelle ligne à la df except: pass SCE_abis_min = dfa_bis[dfa_bis['SCE'] == dfa_bis['SCE'].min()]#on cherche la meilleure valeur de a_bis pour le a donné SCE_abis_min['a'] = a dfa = pd.concat([dfa,SCE_abis_min], ignore_index=True) SCE_a_abis_min = dfa[dfa['SCE'] == dfa['SCE'].min()] #on cherche la combinaison de a et a' qui donne le plus petit SCE SCE_a_abis_min['r'] = r globals()[vent[i]] = pd.concat([globals()[vent[i]],SCE_a_abis_min], ignore_index=True) #matrice indiquant le meilleur a et a' pour une position donnée (une matrice par vitesse de vent) a_best = math.floor(globals()[vent[i]]['a'].iloc[0]) #on prend le meilleur a qui a été trouvé itérations suivantes a_bis_best = math.floor(globals()[vent[i]]['a_bis'].iloc[0]) else : # on va ici se baser sur les valeurs précédentes for z in range (a_best-1,101,1): #on part de la valeur précédente avec une sécurité dfa_bis = pd.DataFrame(columns=['a_bis','CL','CD','Tforce','Qforce','Qmoment','Tmoment','SCE'])#je le fais en deux fois parce que a et a' ne vont pas toujours varier sur la même gamme de valeur ni sur le même nombre de valeurs for z_bis in range (a_bis_best-1,z+1,1): #on part de la valeur précédente avec une sécurité try: a = z/100 #obligé de faire ça comme ça: la boucle fonctionne pas si je veux la faire tourner avec des nombres à virgule a_bis = z_bis/100 phi = math.degrees(math.atan((U*(1-a))/(w*r*(1+a_bis)))) cosphi = math.cos(math.radians(phi)) sinphi = math.sin(math.radians(phi)) #v_res = ((1-a)*U)/(sinphi) v_res = (((1-a)*U)**2)/(w*r*(1+a_bis))**2 alpha = phi - b smallest_gap = 3000 for x in df['alpha']: #boucle pour trouver la valeur du tableau la plus proche du alpha trouvé if abs(alpha - x)<smallest_gap: smallest_gap = abs(alpha - x) alpha_tableau = x CL = df[df['alpha'] == alpha_tableau]['Cl'].iloc[0]#on prend la valeur de CL pour l'alpha correspondant. Le [0] permet d'avoir un chiffre et pas une dataframe. CD = df[df['alpha'] == alpha_tableau]['Cd'].iloc[0] Tforce = 0.5*B*rho*(v_res**2)*(CL*cosphi+CD*sinphi)*c*r Qforce = 0.5*B*rho*(v_res**2)*(CL*sinphi-CD*cosphi)*c*((r**2)/2) Tmoment = 4*np.pi*rho*(U**2)*(1-a)*a*((r**2)/2) Qmoment = 4*np.pi*rho*U*(1-a)*a_bis*w*((r**4)/4) SCE = (((Tforce-Tmoment)**2)+((Qforce-Qmoment)**2))**(1/2) dfa_bis_new = pd.DataFrame(data=np.array([[a_bis,CL,CD,Tforce,Qforce,Qmoment,Tmoment,SCE]]), columns=['a_bis','CL','CD','Tforce','Qforce','Qmoment','Tmoment','SCE']) dfa_bis = pd.concat([dfa_bis,dfa_bis_new], ignore_index=True) #on ajoute une nouvelle ligne à la df except: pass SCE_abis_min = dfa_bis[dfa_bis['SCE'] == dfa_bis['SCE'].min()] SCE_abis_min['a'] = a dfa = pd.concat([dfa,SCE_abis_min], ignore_index=True) SCE_a_abis_min = dfa[dfa['SCE'] == dfa['SCE'].min()] #on cherche la combinaison de a et a' qui donne le plus petit SCE SCE_a_abis_min['r'] = r globals()[vent[i]] = pd.concat([globals()[vent[i]],SCE_a_abis_min], ignore_index=True) #matrice indiquant le meilleur a et a' pour une position donnée (une matrice par vitesse de vent) a_best = math.floor(globals()[vent[i]]['a'].iloc[0]) #on prend le meilleur a qui a été trouvé itérations suivantes a_bis_best = math.floor(globals()[vent[i]]['a_bis'].iloc[0])

plt.plot(huit['r'], huit['a'], color = 'red') plt.plot(huit['r'], huit['a_bis'], color = 'blue') plt.title("Values of a and a' along the blade for a U_wind of 8m/s") plt.xlabel('Distance to the center of the wind turbine [m]') plt.legend(['a','a prim']) plt.show() plt.plot(vingtdeux['r'], vingtdeux['a'], color = 'red') plt.plot(huit['r'], huit['a_bis'], color = 'blue') plt.title("Values of a and a' along the blade for a U_wind of 22m/s") plt.xlabel('Distance to the center of the wind turbine [m]') plt.legend(['a','a prim']) plt.show()

#determination of the total torque for each wind speed dfQ = pd.DataFrame(columns=['Qtot','Wspd']) for i in range (len(vent)): for j in range (len(data_bis['Rayon'])): if j == 0: Qtot = (globals()[vent[i]]['Qmoment'][j])*(globals()[vent[i]]['r'][j]) #value of the total torque approximated for a given section of the blade else: section = globals()[vent[i]]['r'][j] - globals()[vent[i]]['r'][j-1] Qtot = Qtot+(globals()[vent[i]]['Qmoment'][j])*section #sum of the total torque for all the sections of the blades dfQonewind = pd.DataFrame(data=np.array([[Qtot,data['Vitesse du vent'][i]]]), columns=['Qtot','Wspd']) dfQ = pd.concat([dfQ,dfQonewind], ignore_index=True) plt.plot(dfQ['Wspd'],dfQ['Qtot']) plt.title("Values of the torque for differents wind speeds") plt.xlabel('Wind speeds [m/s]') plt.ylabel('torque [Nm]') plt.show() #determination of the power for each wind speed Power = pd.DataFrame(np.zeros((len(vent), 1))) Power['Wspd'] = data['Vitesse du vent'] Power['P'] = dfQ['Qtot']*data['vitesse de rotation'] plt.plot(Power['Wspd'],Power['P']) plt.title("Values of the power produced at different wind speeds") plt.xlabel('Wind speeds [m/s]') plt.ylabel('Electrical power [W]') plt.show()

donnees = pd.read_csv('/work/SWINOUJSCIE-hour-20211124-164632.csv') #The interesting columns are the date, time and Wdir (wind direction) and Wspd (wind speed) donnees_utiles = donnees.iloc[:,[0,1,43,46]] donnees_utiles.head()

Wspd = donnees_utiles['Wspd (m/s)'] vent = Wspd ventbon= sorted(vent) valun = np.ones(8760) DF = pd.DataFrame(valun, columns = ['nombre']) DF['vtrier'] = ventbon for i in range (len(valun)): if i == 0: cumuler=pd.DataFrame(columns=['nombre','vent']) a=1 else: if DF['vtrier'][i]==DF['vtrier'][i-1]: a=a+1 else: dv = pd.DataFrame(np.array([[a,DF['vtrier'][i-1]]]), columns = ['nombre','vent']) cumuler = pd.concat([cumuler,dv], ignore_index=True) a=1 dv = pd.DataFrame(np.array([[a,DF['vtrier'][8759]]]), columns = ['nombre','vent'])#on ajoute la dernière valeur de la boucle cumuler = pd.concat([cumuler,dv], ignore_index=True) cumuler['proba']=cumuler['nombre']/len(valun) cumuler['probacumul'] = np.zeros(len(cumuler)) for i in range(len(cumuler)): cumuler['probacumul'][i]=sum(cumuler['proba'][0:i]) cumuler['X'] = np.log(cumuler['vent']) cumuler['Yf(x)'] = np.log(-np.log(1-cumuler['proba'])) cumuler['YF(x)'] = np.log(-np.log(1-cumuler['probacumul'])) cumul = cumuler.drop((cumuler).index[[0]]) #on se débarrasse des valeurs qui tendent vers l'infini from scipy import stats slope, intercept, r_value, p_value, std_err = stats.linregress(cumul['X'], cumul['YF(x)']) def predict(x): return slope * x + intercept fitLine = predict(cumul['X']) plt.plot(cumul['X'], fitLine, c='r') plt.plot(cumul['X'], cumul['YF(x)']) plt.show() print('The slope is '+str(slope)+' and the horizontal axis is ' + str(intercept)) k = slope print('k = '+ str(k)) c = math.exp(-intercept/k) print('c = '+ str(c)) proba_v_y = np.exp(-((cumuler['vent']/c)**k)) cumuler['heures'] = proba_v_y*365*24 weibully=(k/c)*((cumuler['vent']/c)**(k-1))*np.exp(-((cumuler['vent']/c)**k)) plt.plot(cumuler['vent'],weibully) plt.xlabel('Wind speed [m/s]') plt.ylabel('Probability P(U)') plt.title('Probability of observing a certain wind speed.') plt.show() plt.plot(cumuler['heures'],cumuler['vent']) plt.xlabel('Number of hours per year') plt.ylabel('Wind speed [m/s]') plt.title('Number of hours per year when at least a certain wind speed is observed.') plt.show()

heurebon=pd.DataFrame(columns=['heure']) for i in range (len(cumuler)): if math.fmod(i,10)==0: heure=cumuler['heures'][i] dh = pd.DataFrame(np.array([[heure]]), columns = ['heure']) heurebon = pd.concat([heurebon,dh], ignore_index=True) heurebon.head()

#Determination d'une monotone de puissance # graphe 1 V = np.arange(0,26,1) P = [0, 0, 0, 30, 120, 380, 780, 1290, 1910, 2735, 3740, 4930, 5300, 5200, 5170, 5150, 5140, 5130, 5125, 5120, 5115, 5100, 5000, 5000, 4980, 4955] plt.plot(V, P) plt.xlabel('Wind speed [m/s]') plt.ylabel('Power [kW]') plt.title('Power curve as a function of wind speed') plt.show() # graphe 2 P_V = {'Vitesse du vent' : [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11.0, 12.0], 'Puissance' : [0, 0, 0, 30, 120, 380, 780, 1290, 1910, 2735, 3740, 4930, 5300] } df_puiss = pd.DataFrame(P_V) P = df_puiss['Puissance'] # Power en kW plt.plot(heurebon, P) plt.xlabel('Number of hours per year') plt.ylabel('Power [kW]') plt.title('Power monotone as a function of the number of hours per year') plt.show() vitesse_moy = Wspd.mean() print('The average wind speed = ' + str(vitesse_moy) + ' m/s')

x = heurebon['heure'] y = P # courbePuissance IntegralTrapeze = intg.trapz(x, y) # kWh (The integrals give kWh since the power is expressed in kW) IntegraltrapezeGWh = IntegralTrapeze/10**6 # GWh print('The enrergy is ' + str(IntegralTrapeze) + ' kWh, wich corresponds to '+ str(IntegraltrapezeGWh) + ' GWh')

# load factor nombre_heure = 8760 # number of hours in one year puissance_moy = IntegralTrapeze/nombre_heure print('The average power produce over a year for an average wind speed of 3,91 m/s is ' + str(puissance_moy) + ' kW') nominal_power = P.max() print('The nominal power of the wind turbine for an average wind speed of 3,91 m/s is ' +str(nominal_power) + ' kW') load_factor = puissance_moy/nominal_power*100 print ('The load factor is ', load_factor, '%')