%WORK WITH TURBINES OF 7 METRES
%FLOW DATA
Q=[0:100:12000]; %m^3/s
nu=12.36*10^-4; %dynamic viscosity (kg/ms)
L=18; %m
g=9.81; %m/s^2
rho=1026; %kg/m^3
%TURBINE DATA & CULVERTS
d=7; %m
A=pi*d^2/4; %m^2
N=19; %number of turbines
L=18; %m
%Flow & per turbine
Qt=Q./N; %m^3/s
vc=Qt./A; %m/s
%SECONDARY LOSSES (TRANSITION FROM FREE FLOW TO THE CULVERTS
%could be important for the final amount of power
%REYNOLDS NUMBER
%to calculate the energy losses in the culverts
Re=(rho*d.*vc)./A;
Tablare=[Qt',Re'];
%Once Reynolds is known there are different ways of calculate f (friction
%factor) and therefore, different ways of calculate the energy losses.
if Re<2000
f=64./Re;
elseif 2000<Re<10^5
f=0.3164./Re.^0.25;
end
%Now f is known, the energy losses can be calculated
hr=(8.*f*L.*(Qt.^2))./(g*(pi^2)*d^5);
Tabshr=[Qt',hr'];
%SECONDARY LOSSES WITH THE APARITION OF THE TURBINE??
%SPEED IN THE TURBINE
%We are going to apply Bernouilli to know what's the speed in the turbine
%For that we know that B1=B2+hr
vt=sqrt((vc.^2)-(hr.*2*g));
%PLOT OF THE RESULTS IN ENERGY LOSSES
% figure
% subplot(2,2,[1,2])
% plot(Q,hr)
% title('Losses in the culvert')% xlabel('flow of the barrier (m^3/s)')% ylabel('total losses (wcm)')% subplot(2,2,3)
% plot(Q,Qt,'--')
% title('Comparisson of flows')% xlabel('flow of the barrier (m^3/s)')% ylabel('flow of the turbine (m^3/s)')% subplot(2,2,4)
% plot(Q,vc,Q,vt,'--')
% title('Comparisson of velocities')% xlabel('flow of the barrier (m^3/s)')% ylabel('velocities (m/s)')%
%Once the speeds in the turbine are known, the studio of the velocities in
%the blades can start
%SPEED TRIANGLES
%now we are going to divide the blades in 0.1 width parts
di=2.8; %m inside diameter
dr=[di/2:0.1:d/2]; %m radius diferential
Ar=pi*(d^2-di^2)/4; %m^2 rotor area
vr=vt.*(A/Ar); %m/s velocity in the rotor
nrpm=100; %rpm revolution per minute
n=nrpm/60; %rps revolution per second
%BLADE VELOCITY
du=n*2*pi.*dr;
% figure
% bar(du,'blue')
cm=vr;
for q= 100
%we are going to work with a flow speed (1m/s) and
%start with the position of the blade. Once we know how the
%blade is gonna be in that flow, we only need to evaluate that
%blade in the other flows.
z=0;
while z<length(dr)
z=z+1;
%VELOCITIES
% beta=atand(wm1/wuf);
% alfa=atand(wm1/cuf);
%ENTRY
alfa1=[2:1:89]; %º all angles posibilities
c1=cm(q)./sind(alfa1); %m/s
cu1=c1.*cosd(alfa1); %m/s
cu1s(z,:)=cu1;
wu1=du(z)-cu1; %m/s (problemas con las dimensiones)
wu1s(z,:)=wu1; %m/s
wm1=cm(q); %m/s
w1=sqrt(wu1.^2+wm1.^2); %m/s
beta1=atand(wm1./wu1); %º
beta1s(z,:)=beta1;
w1s(z,:)=w1; %m/s all relative velocities at the entry
T1=[alfa1',beta1',w1']; %º º m/s
%SELECTION OF THE FINAL ANGLES (ENTRY)
bs=beta1s;
bs(bs<0)=99.99;
entryangles=min(bs');
[x,y]=find(bs'==entryangles(z));
good_alfa=alfa1(x);
entryalfas(z)=good_alfa;
%EXIT
alfa2=[2:1:89]; %º
c2=cm(q)./sind(alfa2); %m/s
cu2=c2.*cosd(alfa2); %m/s
wu2=du(z)-cu2; %m/s (supongo que aquí también)
wm2=cm(q); %m/s
w2=sqrt(wu2.^2+wm2.^2); %m/s
beta2=atand(wm2./wu2); %º
beta2s(z,:)=beta2;
w2s(z,:)=w2; %m/s
T2=[alfa2',beta2',w2']; %º º m/s
%SELECTION OF THE FINAL ANGLES (EXIT)
bes=beta2s;
exitangles=max(bes');
[x,y]=find(bs'==exitangles(z));
exit_alfa=alfa2(x);
exitalfas(z)=exit_alfa;
%NEW TRIAL WITH THE ANGLES
finalangles=[entryangles,exitangles];
new_alfa1=entryalfas;
new_c1=cm(q)./sind(new_alfa1);
new_cu1=new_c1.*cosd(new_alfa1);
%new_cu1s(z,:)=new_cu1;
new_wu1=du(z)-new_cu1(z);
new_wu1s(z)=new_wu1;
new_wm1=cm(q);
new_w1=sqrt(new_wu1.^2+new_wm1.^2);
new_w1s(z)=new_w1;
new_alfa2=exitalfas;
new_c2=cm(q)./sind(new_alfa2);
new_cu2=new_c2.*cosd(new_alfa2);
%new_cu2s(z,:)=new_cu2;
new_wu2=du(z)-new_cu2(z);
new_wu2s(z)=new_wu2;
new_wm2=cm(q);
new_w2=sqrt(new_wu2.^2+new_wm2.^2);
new_w2s(z)=new_w2;
new_ws=[new_w1s',new_w2s'];
new_ws2=[new_w1s'.^2,new_w2s'.^2];
deltaP=(rho.*((new_w1s(z)^2)-(new_w2s(z)^2)))./2;
deltaPs(z)=deltaP;
deltamax=(6.5*g*rho)+101325;
pressure=deltamax-deltaPs;
%m/s
while pressure(z)<10^3
entryalfas(z)=entryalfas(z)-1
new_alfa1=entryalfas;
new_c1=cm(q)./sind(new_alfa1);
new_cu1=new_c1.*cosd(new_alfa1);
%new_cu1s(z,:)=new_cu1;
new_wu1=du(z)-new_cu1(z);
new_wu1s(z)=new_wu1;
new_wm1=cm(q);
new_w1=sqrt(new_wu1.^2+new_wm1.^2);
new_w1s(z)=new_w1;
new_alfa2=exitalfas;
new_c2=cm(q)./sind(new_alfa2);
new_cu2=new_c2.*cosd(new_alfa2);
%new_cu2s(z,:)=new_cu2;
new_wu2=du(z)-new_cu2(z);
new_wu2s(z)=new_wu2;
new_wm2=cm(q);
new_w2=sqrt(new_wu2.^2+new_wm2.^2);
new_w2s(z)=new_w2;
deltaP=rho.*((new_w1s.^2)-(new_w2s.^2))./2
deltamax=(6.5*g*rho)+101325; %min total depth of 14 metres
pressure(z)=deltaP(z)-deltamax;
if pressure(z)>10^3
cuf=(new_cu2(z)-new_cu1(z))/2;
end
end
end
end
