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
z=2;                                                                                        %Number of blades
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 
w=2*U*TSR/D;                                                                                %rotation rad/s
[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);
            phirad(i)=atan(((1-aa(i))*U)/((1+ar(i))*w*r(i)))  ;                                %Phi computation with a(i) and a'(i)gh [rad]
            phi(i)=phirad(i)*180/pi;
            Tw(i)=phi(i)-alpha(i);                                                       
            phi(i)=phirad(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.

Delete this entry (admin only).