%% Turbine Design Parameters
R = 1; % Turbine Radius
N = 1; % Number of Blades
c = 1; % Chord Length
L = 1; % blade Height
%% Operating Parameters
Vo = 7; % Free Stream Velocity (m/s)
n=36; % Number of Streamtubes ( 180 / 5 = 36 ) ???
w = 8; % Angular Velocity (rad/s)
Xt = w*R/Vo; %Tip speed ratio
Ao = 3; % initial angle of attack (rad)
%% CONSTANTS FOR STANDARD AIR CONDITIONS AND SEA LEVEL
kv = 1.4607*(exp(-5)); % Kinematic viscosity at 15oC [m^2/s] ???
rho = 1.225; % air density (Standard density at sea level)[kg/m^3] at sea level ????
%% Process
S = 2*L*R; %Swept area [m^2]
theta = -90;
theta2 = -89;
st = 1;
for st = 1:1:36
if theta < 0
theta = theta + 360;
end
if theta2 < 0
theta2 = theta2 + 360;
end
if theta - 360 < 0 % ?????
theta = theta - 360 + (179 / 36) ;
else
theta = theta + (179 / 36);
end
if theta2 - 360 < 0
theta2 = theta2 - 360 - (179 / 36);
else
theta2 = theta2 - (179 / 36);
end
end
thetau(st) = theta * pi / 180;
thetad(st) = (theta2) * pi / 180;
Vu = Vo;
%%---------------UPSTREAM CALCULATION---------------------
i = 0; %initialize the counter value
while (i~=n)
i = i+1;
% Wu=local resultant air velocity
Wu = sqrt( Vu^(2)*( (Xt - sin ( thetau(i)))^(2) + (cos ( thetau(i)))^(2)));
Reb = Wu*c/kv; %Local Blade Reynolds number
%% Airfoil
% Type of airfoil
typeNACA = '2412';
% Extracrt values from type of airfoil string
Minit = str2double(typeNACA(1));
Pinit = str2double(typeNACA(2));
Tinit = str2double(typeNACA(3:4));
% Actual percentage values of airfoil properties
M = Minit/100;
P = Pinit/100;
T = Tinit/100;
%%
%if NACA == 15
%[A1, Cl1, Cd1] = NACAfinder15(Reb);
%else
%[A1, Cl1, Cd1] = NACAfinder21(Reb);
%end
A1=0:5:180;
cd1=1:1:36;
cl1=0.1:0.03:0.7;
costh = cos(thetau(i));
cosao = cos(Ao);
sinth = sin(thetau(i));
sinao = sin(Ao);
A =asin(costh*cosao-(Xt-sinth)*sinao/sqrt((Xt-sinth)^2+(costh^2)));
Cd = interp1(A1,cd1,A,'linear');
Cl = interp1(A1,cd1,A,'linear');
switch A
case 0<A && A<60
Cd = 0.5;
Cl = 0.5;
case 60<A && A<140
Cd = (200 - 2 * Ao) / 80;
Cl = (200 - 2 * Ao) / 80;
case 140<A && A<160
Cd = (A - 173.33) / 26.66;
Cl = (A - 173.33) / 26.66;
case 160<A && A<180
Cd = (A - 180) / 80;
Cl = (A - 180) / 80;
end
au = 1.01; % velocity induction factor upstream
newau = 1; % initialize, au must be different from newau ??????
while (au-newau)>(10^(-3))
Vu = Vo*(au); %Vu =velocity upstream of wind turb
X = R*w/Vu; %Local Tip speed ratio
if ((A) < 0)
A = -1 * A;
Cl = -1 * Cl;
end
% Wu=local resultant air velocity
Wu = sqrt( Vu^(2)*( (X - sin ( thetau(i)))^(2) + (cos ( thetau(i)))^(2)));
% interp1 is 1-D data interpolation for example x goes 1 to 5 , y has 6
% point and also you can define 18 points on that figure, those 6 points
% are main points for your graph but also the other 18 points are wanted to
% show on that figure.
%% Continue
% Cn = normal coefficient, % Ct = tangential coefficient
Cnu = Cl*cosd (A) + Cd*sind (A); %note that A is in degrees
Ctu = Cl*sind (A) - Cd*cosd (A);
% fup = function to find interference factor (au)
g=@(thetau) (abs(sec (thetau)).*(Cnu.*cos(thetau) - thetau) - Ctu.*sin(thetau)).*(Wu./Vu).^2;
for g = 0:10
fup(g) = ((Abs(1 / ((Cos((-pi / 2) + (g * pi /10)))))) * (cnu * Cos((-pi / 2) + (g * pi / 10)) - ctu * Sin((-pi / 2)+ (g * pi / 10))));
actfup = (N * c* y / (8 * pi * R ))* ((pi / 30) * ((fup(0) + fup(10)) + (4 * (fup(1) + fup(3) + fup(5) + fup(7) + fup(9))) + 2 * (fup(2) + fup(4) + fup(6) + fup(8))));
end
newau = pi/(actfup+pi); % new interference factor value for the next iteration process
Fnu = (c*L/S)*Cnu*(Wu/Vo)^2; % normal force coefficient
Ftu = (c*L/S)*Ctu*(Wu/Vo)^2; % tangential force coefficient
Tup(i) = 0.5*rho*c*R*L*Ctu*(Wu^2);
av_tup_i = (1 / 3) * ((Tup(1) + Tup(36)) + 4 * (Tup(3) + Tup(5) + Tup(7) + Tup(9) + Tup(11) + Tup(13) + Tup(15) + Tup(17) + Tup(19) + Tup(21) + Tup(23) + Tup(25) + Tup(27) + Tup(29) + Tup(31) + Tup(33) + Tup(35)) + 2 * (Tup(2) + Tup(4) + Tup(6) + Tup(8) + Tup(10) + Tup(12) + Tup(14) + Tup(16) + Tup(18) + Tup(20) + Tup(22) + Tup(24) + Tup(26) + Tup(28) + Tup(30) + Tup(32) + Tup(34)));
av_Tup = N*(av_tup_i)/(2*pi); % up stream average torque
av_Cqu = av_Tup/(0.5*rho*S*R*Vo^2);
Cpu = av_Cqu*Xt; % Upstream power coefficient
Auvector (i) = A; % Store angle of attack value
auvector (i) = newau; % Store au value in a vector
end
%%---------------DOWNSTREAM CALCULATION---------------------
j = n+1;
flag =0;
i = 0; %initialize the counter
while (j~=1)
j = j-1;
i = i+1; % interference factor downstream.
ad = 1.01; % velocity induction factor upstream
newad = auvector(j); % initialize, au must be different from newau
while (ad-newad)> 0.001 %Iterative process to find ad
ad = newad;
Ve = Vo*((2*auvector(j))-1); %Ve = air velocity inside cylinder
Vd = Ve*ad;
if vd > 0
X = r * w / vd;
else
X = ve * 1.01;
end
Wd = sqrt ( Vd^2*( (X -sin (thetad(i)))^2 + (cos (thetad(i)))^2));
Reb = Wd*c/kv;
if NACA == 15
[A1, Cl1, Cd1] = NACAfinder15(Reb);
else
[A1, Cl1, Cd1] = NACAfinder21(Reb);
end
costh = cos(thetad(i));
cosao = cos(Ao);
sinth = sin(thetad(i));
sinao = sin(Ao);
A2 = asin((costh*cosao -(X-sinth)*sinao)/sqrt((X-sinth)^2+(costh^2)));
neg = 0;
if ((A2) < 0)
neg = 1
end
A2 = abs(A2*180/pi);
Cl = interp1 (A1, Cl1, A, 'linear');
Cd = interp1 (A1, Cd1, A, 'linear');
switch A2
case 0<A && A<60
cd = 0.5;
Cl = 0.5;
case 60<A && A<140
cd = (200 - 2 * Ao) / 80;
Cl = (200 - 2 * Ao) / 80;
case 140<A && A<160
cd = (Ao - 173.33) / 26.66;
Cl = (Ao - 173.33) / 26.66;
case 160<A && A<180
cd = (Ao - 180) / 80;
Cl = (Ao - 180) / 80;
end
if (neg==1)
A2 = -1*A2;
Cl =-1*Cl;
Cd =1*Cd;
end
Cnd = Cl*cosd (A2) + Cd*sind (A2);
Ctd = Cl*sind (A2) - Cd*cosd (A2);
g=@(thetad) (abs(sec (thetad)).*(Cnd.*cos(thetad) - Ctd.*sin(thetad)).*(Wd./Vd).^2);
for g1 = 0 : 10
y = integral (@g, 91*pi/180, 269*pi/180);
fdw(g1) = ((abs(1 / (cos((-pi / 2) + (g1 * pi / 10))) ^ (-1))) * (Cnu * cos((-pi / 2) + (g1 * pi / 10)) - Ctu * sin((-pi / 2) + (g1 * pi / 10))));
actfdw = (N * c * y / (8 * pi * r)) * ((pi / 30) * ((fdw(0) + fdw(10)) + (4 * (fdw(1) + fdw(3) + fdw(5) + fdw(7)+ fdw(9))) + 2 * (fdw(2) + fdw(4) + fdw(6) + fdw(8))));
end
if (flag ==0)
newad = pi/(fdw+pi);
end
if (newad<0.01)
warning('newad< 0.01 at theta = %d and A = %d', (thetad(i)*180/pi), A);
if (i>1)
newad = advector(i-1);
else
newad = auvector (i);
end
flag = 1;
end
end
Advector (i) = A;
advector (i) = newad;
Fnd (i) = (c*L/S)*Cnd*(Wd/Vo)^2; % normal force coefficient
Ftd (i) = (c*L/S)*Ctd*(Wd/Vo)^2; % tangential force coefficient
Tdw (i) = 0.5*rho*c*R*L*Ctd*Wd^2;
end
end
% Average upstream torque
ts4 = f_trapezoidal_integration_s (thetad, Tdw);
av_Tdw = N*(ts4)/(2*pi); % upstream average torque age torque
% Average torque coefficient
av_Cqd = av_Tdw/(0.5*rho*S*R*Vo^2);
Cpd = av_Cqd*Xt; % downstream power coefficient
Cpt = Cpd+Cpu; % Total power coefficient
av_T = av_Tup + av_Tdw; % Total average torque [Nm]
