# Demo entry 6662308

MATLAB AIR LIFT

Submitted by RIC on Nov 25, 2017 at 15:39
Language: Matlab. Code size: 6.6 kB.

```clc;
clear all ;
%% input parametres
D=0.9;                                                                                      %Diameter
U=10;                                                                                       %Wind velocity
TSR=7.5;                                                                                    %Tip speed ratio
ro=1.2;                                                                                     %Wind density
Cpd=0.5;                                                                                    %Set up required value for Cpd
nel=10;                                                                                     %number of blade of elements

%% Preparation
r=zeros(1,nel);                                                                             %set up r
for i = 1:nel                                                                               %radius of each element
r(i)=D/(2*nel)*i-D/(4*nel);
end
dr=D/(2*nel);                                                                               %Width of each element
[Cl,Cd,alpha]=lift_drag_coefNACA                                                                %optimal values for Cl,Cd and alpha
Lc=[0.1 0.063 0.043 0.032 0.026 0.021 0.019 0.017 0.015 0.014];                             %First Excesice
%% Setting vectors
A=(pi*D^2)/4;                                                                               %Area of rotor                                                                           %number of blade of elements
aa=0.33*ones(1,nel);                                                                        %Set up a (first guess)
ar=zeros(1,nel);                                                                            %Set up a' (first guess)
aa0=zeros(1,nel);                                                                           %Set up a(i-1) firstly to start iteration
ar0=zeros(1,nel);                                                                           %Set up a'(i-1) firstly to start iteration
phi=zeros(1,nel);                                                                           %Set up vectors for calculation
Ca=zeros(1,nel);
Cr=zeros(1,nel);
sigm=zeros(1,nel);
Bep=zeros(1,nel);
dT=zeros(1,nel);
dM=zeros(1,nel);
alpha=alpha*ones(1,nel);
Cl=Cl*ones(1,nel);
Cd=Cd*ones(1,nel);
M=0;
for i=1:nel                                                                                 %Compute sigma for each elemment
sigm(i)=z*Lc(i)/(2*pi*r(i));
end
Cp=0;                                                                                       %Set up value for Cp

i=1;                                                                                        %New i starting value for loop
%% Main calculations

for i=1:nel                                                                            %For each blade element
while abs(aa0(i)-aa(i))>10^(-10) || abs(ar0(i)-ar(i))>10^(-10)                          %Condition for iteration
aa0(i)=aa(i);
ar0(i)=ar(i);
Tw(i)=phi(i)-alpha(i);
if alpha(i)<-180
break
end
if alpha(i)>180
break
end
Ca(i)=Cl(i)*cos(phi(i))+Cd(i)*sin(phi(i));                                      %Ca computation
Cr(i)=Cl(i)*sin(phi(i))-Cd(i)*cos(phi(i));                                      %Cr computation

F=2/pi*acos(exp(-z*((D/2)-r(i))/(2*r(i)*sin(phi(i)))));
if aa0(i)<0.2                                                                   %Condition for guessed a(i)
aa(i)=((4*F*sin(phi(i))*sin(phi(i))/(Ca(i)*sigm(i)))+1)^(-1);                          %Computed new values for a(i+1)
ar(i)=((4*F*sin(phi(i))*cos(phi(i))/(Cr(i)*sigm(i)))-1)^(-1);                 %Computed new values for a'(i+1)
else

K=4*F*sin(phi(i))*sin(phi(i))/(Ca(i)*sigm(i));

if ((K*(1-2*0.2)+2)^2+4*(K*0.2^2-1))<0
break
end
aa(i)=0.5*(2+K*(1-2*0.2)-sqrt((K*(1-2*0.2)+2)^2+4*(K*0.2^2-1)));
ar(i)=((4*F*sin(phi(i))*cos(phi(i))/(Cr(i)*sigm(i)))-1)^(-1);

end

%New value of a'(i+1) is raised by delta and becomes a'(i) for next loop

end

Bep(i)=aa(i)/(1+ar(i))*4*sin(phi(i));                                                 %Computation of Bep
Lc(i)=Bep(i)*2*pi*U/(z*Cl(i)*w);                                                      %New values of Cord length
dT(i)=F*z*Ca(i)*0.5*ro*Lc(i)*U^2*(1-aa(i))^2*dr/(sin(phi(i)).^2);                       %Torque of each element
dM(i)=F*z*Cr(i)*0.5*ro*Lc(i)*U*(1-aa(i))*r(i)*w*(1+ar(i))*dr/(sin(phi(i))*cos(phi(i))); %Momentum of each element
M=M+dM(i)*r(i);                                                                       %raise of Momentum for whole turbine
if i<nel
aa(i+1)=aa(i);
ar(i+1)=ar(i);
end
end
T=sum(dT)
P=M*w;
Cp=P/(0.5*ro*U^3*A);
%% Output

%1. Maximum power coefficient
Cp
%2. Maximum power output, P [W]
P
%3. Thrust force, T [N]
T
%4. Plot showing the torque and thrust dirstribution versus radius
figure(1)
plot(r/D/2,dT,r/D/2,dM)
xlabel('r/R')
legend('Thrust','Torque')
%5. Plot showing the cord length versus radius
figure(2)
plot(r/D/2,Lc)
xlabel('r/R')
legend('Chord length','Location','northeastoutside');
%6. Plot showing the Cl/Cd ratio for each of the blade elements
figure(3)
plot(r/D/2,Cl./Cd)
xlabel('r/R')
legend('Cl/Cd')

%6. Plot showing the angle of attack for each of the blade elements
figure(5)
plot(r/D/2,alpha)
xlabel('r/R')
legend('alpha')

```

This snippet took 0.03 seconds to highlight.

Back to the Entry List or Home.