%% Exercise sur le calcul du facteur d'induction d’une turbine éolienne

clear all;
close all;
clc;

% Une turbine éolienne à axe horizontal est exploitée avec les conditions suivantes
% vitesse moyenne du vent :U_{\infty } =14\;\frac{m}{s}
% densité de l’air : \rho =1\ldotp 205\;\frac{\textrm{kg}}{m^3 }
% coefficient de puissance : C_P =0\ldotp 29
% puissance nominale du générateur électrique : P=2\ldotp 2\textrm{MW}

U_inf=14 %[m/s]
rho=1.205 %[kg/m^3]
Cp=0.29
P=2.2E6 %[W]

%% Questions
% Calculer le radius de la turbine que, en correspondance de la valeur de vitesse moyenne du vent, permet le couplage avec le générateur électrique;
% Calculer la puissance maximale que la turbine peut produire (en correspondance de la valeur de vitesse moyenne du vent) pour le C_P correspondent à la limite de Lanchester-Betz;
% Déterminér la valeur de la vitesse du vent en correspondance du disque actuateur de la turbine U_D.
%
%% Solution
% Question 1
% A travers la définition du coefficient de puissance, C_P =\frac{P}{\frac{1}{2}\rho \;U_{\infty }^3 A_D }, il est possible de déterminer le diamètre de la turbine.
% A_D =\frac{\pi D^2 }{4}=\frac{P}{\frac{1}{2}\rho U_{\infty }^3 A_D }
% D=2\sqrt{\frac{P}{\pi \frac{1}{2}\rho U_{\infty }^3 C_P \;}\;}

D=2*sqrt(P/(0.5*pi*rho*(U_inf^3)*Cp)) %[m]

%
% Question 2
% La valeur maximale de la limite de Lanchester-Betz est, pour une turbine à axe horizontal, est 16/27, donc la puissance maximale est  simplement
% P_{\mathrm{max}} =\frac{1}{2}C_P^{\mathrm{max}} \rho U_{\infty }^3 A_D

A_D=pi*D^2/4; %[m^2]
Cp_max=16/27;
P_max=0.5*Cp_max*rho*U_inf^3*A_D %[W]

% Question 3
% La vitesse en correspondance du disque actuateur, U_D, peut être déterminée à travers la définition du facteur d’induction a:
% U_D =U_{\infty } \left(1-a\right)
% Et la valeur du facteur d’induction peut être déterminée à travers la relation entre C_P \;et a:
% C_P =4{a\left(1-a\right)}^2
% Donc, on peut determiner les racines de la fonction C_P -4{a\left(1-a\right)}^2 =0 pour le C_P =0\ldotp 29

a0=0; %valeur initiale pour la solution de l'equation
a = fsolve(@(a) Cp-4*a*(1-a)^2,a0)

% Donc, la vitesse du vent en correspondence du disque est:
% U_D =U_{\infty } \left(1-a\right)

U_D=U_inf*(1-a) %[m/s]