%% Exercise sur la détermination des courbes couple-vitesse (C-Ω) et couple puissance (P-Ω) d’une turbine éolienne

clear all;
close all;
clc;

% Une turbine éolienne à axe horizontal a un rayon R=38\ldotp 2m

R=38.2; %[m]
Ad=pi*R^2; %[m^s] area Ad de la turbine

% En connaissant la courbe Cp-\lambda nous souhaitons déterminer la caractéristique couple-vitesse de rotation de l’arbre de la turbine (C-Ω) dans une plage de variation de vitesse du vent  3\frac{m}{s}\le U_{\infty } \le 12\frac{m}{s}.

Uinf=3:1:12; %[m/s]
b=length(Uinf);

% La densité de l’air est \rho =1\ldotp 205\;\frac{\textrm{kg}}{m^3 }

rho=1.205 %[kg/m^3]

% La discrétisation de la courbe Cp-Ω  est la suivante

%Cp-Lamda discretization
lambda=[1 2 3 3.5 3.75 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 10 11 12 13 14 15];
Cp=[0.015 0.06 0.16 0.27 0.33 0.379 0.437 0.483 0.506 0.523 0.537 0.541 0.538 0.535 0.532 0.527 0.515 0.49 0.449 0.397 0.322 0.242];
a=length(lambda);

plot(lambda,Cp);
xlabel('\lambda')
ylabel('Cp')
title('Courbe CP-\lambda')

%% Questions
% Calculer la puissance et le couple de la turbine en correspondance une vitesse du vent U_{\infty } =5\frac{m}{s} et \lambda =8
% Déterminer les courbes couple-vitesse (C-\omega) et couple puissance (P-\omega) en correspondance de la plage de variation de vitesse du vent  3\frac{m}{s}\le U_{\infty } \le 12\frac{m}{s}.
%
%% Solution
% Pour chaque valeur de vitesse du vent U_{\infty } et chaque valeur de \lambda, il est possible d’évaluer la valeur correspondante de la vitesse de rotation de l’arbre de la turbine \omega par l’équation du "tip-speed ratio":
% \omega =\frac{\lambda U_{\infty } }{R}

omega=(8*5)/R %[rad/sec]

% Avec l'expression de la puissance de la turbine, et la conneissance du  C_P =f\left(\lambda \right), on a:
% P=\frac{1}{2}\rho \;A_D {U_{\infty } }^3 C_P

Pout=0.5*rho*Ad*(5^3)*Cp(13) %[W]

% Le couple correspondant est calculé à l’aide de l’équation suivante
% C=\frac{P}{\omega }

Couple=Pout/omega   %[Nm]

% Changement de variable \omega =f\left(\lambda ,U_{\infty } \right) avec l'équation du "tip-speed ratio"

for i=1:b
    for j=1:a
        omega(i,j)=lambda(j)*Uinf(i)/R;
    end
end

% Calcul de la puissance mechanique de sortie de la turbine pour les valeurs de U_{\infty } et \lambda

for i=1:b
    for j=1:a
        Pout(i,j)=0.5*rho*Ad*Cp(j)*(Uinf(i)^3);
    end
end

% Calcul du couple de la turbine pour les valeurs de U_{\infty } et \lambda

for i=1:b
    for j=1:a
        Couple(i,j)=Pout(i,j)/omega(i,j);
    end
end

% Plot des caracteristiques C-\omega

figure(2)
for i=1:b
    plot(omega(i,:),Couple(i,:));
    hold on
end
xlabel('\omega [rad/sec]')
ylabel('Couple [Nm]')
title('Courbe Couple-\omega')

% Plot des caracteristiques P-\omega

figure(3)
for i=1:b
    plot(omega(i,:),Pout(i,:));
    hold on
end
xlabel('\omega [rad/sec]')
ylabel('Puissance [W]')
title('Courbe Puissance-\omega')