clear all
clc
rpms = [1600, 1800, 2000, 2200, 2400, 2600, 2800, 3000];
dz = 0.2; % Diferencia de cotas entre la entrada y la salida = 20 cm
rho = 1e3; % Densidad = 1000 kg/m3
g = 9.8; % Aceleracio´n de la gravedad
D = 0.2; % Longitud caracter´ıstica del impulsor
% Creamos vectores vac´ıos para almacenar los datos:
vecq = []; vech = []; vecwu = []; vecwt = []; veceta = [];
% CAUDAL ALTURA POT UTIL POT EJE RENDIMIENTO
for i=length(rpms):-1:1
% Definimos el nombre de los ficheros que hay que leer
run = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1162238/robainutile.txt',...
'DecimalSeparator',',') %##### MODIFIED TO RUN HERE
% Convertimos de RPM a RAD/S
run(:,1) = run(:,1)*2*pi/60;
% Convertimos de BAR a PASCAL
run(:,2:3) = run(:,2:3)*1e5;
% La columna 4, representa el par en SI => no se toca
% Convertimos de L/MIN a M3/S
run(:,5) = run(:,5)*1e-3/60;
% Definimos un vector 25 caudal con 500 puntos:
q = linspace(0,max(run(:,5)),500);
% Ordenamos los datos en orden Q creciente
[aux, ind] = sort(run(:,5));
run = run(ind,:);
% Calculamos la altura manome´trica:
gHm = (run(:,3)-run(:,2))/rho + g*dz;
% Ajuste de gHm con una funcion polinomica de segundo orden
[p,S,mu] = polyfit(run(:,5),gHm,2);
% Interpolamos gHm para valores de caudal q:
[y,delta] = polyval(p,q',S,mu);
% Almacenamos los resultados en un vector vecq y vech
vecq = [vecq ; q']; vech = [vech; y];
% Calculamos la potencia útil:
Wu = rho*run(:,5).*gHm;
% Ajustamos y la interpolamos
[pwu,Swu,muwu] = polyfit(run(:,5),Wu,3);
Wufit=polyval(pwu,q',Swu,muwu);
%Almacenamos en vector
vecwu = [vecwu ; Wufit];
% Calculamos la potencia en el eje
Wt = run(:,1).*run(:,4);
% y la interpolamos
[pwt,Swt,muwt] = polyfit(run(:,5),Wt,1);
Wtfit=polyval(pwt,q',Swt,muwt);
% Almacenamos los resultados en un vector vecwt
vecwt = [vecwt ; Wtfit];
% Calculamos rendimiento
eta = Wu./Wt*100;%rho*q.*y' ./Wtfit'*100;
% Interpolamos rendimiento
[peta,Seta,mueta] = polyfit(run(:,5),eta,2);
etafit=polyval(peta,q',Seta,mueta);
% Almacenamos los resultados en un vector veceta
veceta = [veceta ; etafit];
% Hacemos las curvas dimensionales
figure(1); hold on;
xlabel('Q (m^3/s)');ylabel('gHm (m^2/s^2)');
plot(run(:,5),gHm,'ob'); hold on;
plot(q,y,'r-');
grid on;box on;
figure(2); hold on;
xlabel('Q (m^3/s)');ylabel('Wu (W)');
plot(run(:,5),Wu,'ob'); hold on;
plot(q,Wufit,'r-');
grid on;box on;
figure(3); hold on;
xlabel('Q (m^3/s)');ylabel('Wt (W)');
plot(run(:,5),Wt,'ob'); hold on;
plot(q,Wtfit,'r-');
grid on; box on;
figure(4); hold on;
xlabel('Q (m^3/s)');ylabel('\eta (%)');
plot(run(:,5),eta,'ob'); hold on;
plot(q,etafit,'r-');
xlabel('Q (m^3/s)');ylabel('\eta (%)');
grid on; box on;
end
m = 101;
n = 100;
ejeQ = linspace(0,max(vecq),m);
ejeWu = linspace(0,max(vecwu),n);
ejeWt = linspace(0,max(vecwt),n);
ejeHm = linspace(min(vech),max(vech),n);
[QQ,HH] = meshgrid(ejeQ,ejeHm);
[QQ,WW] = meshgrid(ejeQ,ejeWt);
ETAqh = griddata(vecq,vech,veceta,QQ,HH,'linear');
ETAqw = griddata(vecq,vecwt,veceta,QQ,WW,'linear');
indETAqh=isnan(ETAqh);
indETAqw=isnan(ETAqw);
figure;
plot3(vecq,vech/g,veceta,'k.')
hold on
ETAqh2=gridfit(vecq,vech,veceta,ejeQ,ejeHm);
ETAqh2(find(indETAqh==1)) = NaN;
surf(QQ,HH/g,ETAqh2);colormap(jet);colorbar;
camlight right; lighting phong; shading interp;
xlabel('Q (m^3/s)'); ylabel('H_m (m)'); zlabel('\eta (%)');
title('Rendimientos usando GRIDFIT')
box on
figure;
plot3(vecq,vecwt,veceta,'k.')
hold on
ETAqw2=gridfit(vecq,vecwt,veceta,ejeQ,ejeWt);
ETAqw2(find(indETAqw==1)) = NaN;
surf(QQ,WW,ETAqw2);colormap(jet);colorbar;
camlight right; lighting phong; shading interp;
xlabel('Q (m^3/s)'); ylabel('W_t (W)'); zlabel('\eta (%)');
title('Rendimientos usando GRIDFIT')
box on
figure;
[C,h] = contourf(QQ,HH/g,ETAqh2,[0 5 10 15 20 25 30 32.5 35]);colormap(jet);colorbar;
text_handle = clabel(C,h,'FontSize',14,'Color','k','Rotation',0,'LabelSpacing',500);
set(text_handle,'BackgroundColor',[1 1 1])
xlabel('Q (m^3/s)'); ylabel('H_m (m)');
title('Isorendimientos usando GRIDFIT')
figure;
[C,h] = contourf(QQ,WW,ETAqw2,[0 5 10 15 20 25 30 32.5 35]);colormap(jet);colorbar;
text_handle = clabel(C,h,'FontSize',14,'Color','k','Rotation',0,'LabelSpacing',500);
set(text_handle,'BackgroundColor',[1 1 1])
xlabel('Q (m^3/s)'); ylabel('W_t (W)');
title('Isorendimientos usando GRIDFIT')
