
% Método de las características
% Este programa no considera la aparición de cavitación

clc
clear all
close all
 
% Datos
Lgp = 3*154.8 ;     % longitud de la galería de presión(m)
Ltf = 1197.25-Lgp ; % longitud de la tubería forzada(m)
Ltdif = 20 ;        % longitud del tubo difusor

L= Lgp + Ltf ;

landa = 0.02 ;      % coeficiente de fricción del conducto
D = 4.4 ;           % diámetro del conducto (m)
% v_0 = 1 ;           % m/s
c_0 = 1300 ;        % m/s

N = 50 ;
inc_t = 1.0/N ;     % Incremento de tiempo
tfin = 50.0 ;
% tfin = 0.1 ;
t_p = (0:1.0/N:tfin);
lt=length(t_p);     % = tfin*N + 1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % nuevos parámetros para imponer el cierre no instantáneo
t_c = 1;           % tiempo de cierre de la válcula
t_iv = L/c_0;       %tiempo de ida y vuelta de la onda
A_e = (pi*D^2)/4;   % Área de entrada de la válvula
A_s = A_e ;         % Área de salida de la válvula
% a_0 = A_s/A_e ;     % Constante adimensionales para t=0
a_0 = 0.01129200682 ; % Para que salga finalmente v_0 próximo a 1
H_g = 400 ;
g = 9.81 ;

v_0 = sqrt((2*g*H_g)/(1/(a_0^2)+landa*L/D));

 n = 0.5;            % (n<1/2) Cierre muy rápido, siempre aparece compresibilidad
% n = 1.0;              % Ley de cierre lineal
% n = 2.0;              % Cierre lento

a = a_0*((t_c-t_p)/t_c).^n ;   % Ley de cierre

% v_0 = sqrt((2*g*H_g)/(1/(a(1)^2)+landa*L/D));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x_p = (0:1.0/N:1);  
lx=length(x_p);     % = N + 1

% Condiciones de iniciales para t=0 (primera fila) y para toda x (todas las columnas)

% p_p(1,1:lx) = 0 ;
% p_p = rand(1,lx);
for i = 1:lx
    p_p(1,i) = -(landa*L/(2*D))*(v_0/c_0)*i/lx ;
end
v_p(1,1:lx) = 1 ;

% A continuación, se resuelven el resto de los puntos

A = [1 1 ; 1 -1];


ca = [];
cb = [];
cc = [];


for i = 2:lt
    j=1;
    p_p(i,j) = 0 ;
    v_p(i,j) = v_p(i-1,j+1) - p_p(i-1,j+1) - landa*L/(2*D)*(v_0/c_0)*inc_t*v_p(i-1,j+1)*abs(v_p(i-1,j+1)) + p_p(i,j) ;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    j=lx;
    if i <= t_c*N
        coefa = (v_0/(2*c_0*(a(i)^2)));
        coefb = 1;
        coefc = -(v_p(i-1,j-1) + p_p(i-1,j-1) + (g*H_g)/(v_0*c_0) - (landa*L*v_0*inc_t)/(2*D*c_0)*(v_p(i-1,j-1))^2);
        
        ca = [ca coefa];
        cb = [cb coefb];
        cc = [cc coefc];
        
        v_p(i,j) = (-coefb + sqrt(coefb^2-4*coefa*coefc))/(2*coefa);
%           v_p(i,j) = (-1 + sqrt(1-4*(v_0/(2*c_0))*(1/(a(i)^2))*(v_p(i-1,j-1) + p_p(i-1,j-1) - landa*L/(2*D)*(-1)*(v_0/c_0)*inc_t*v_p(i-1,j-1)*abs(v_p(i-1,j-1)) + g*H_g/(v_0*c_0))))/(2*(-v_0/(2*c_0))*(1/(a(i)^2))) ;
%           v_p(i,j) = (-1 + sqrt(1-4*(-v_0/(2*c_0))*(1/(a(i)^2))*(v_p(i-1,j-1) + p_p(i-1,j-1) - landa*L/(2*D)*(v_0/c_0)*inc_t*v_p(i-1,j-1)*abs(v_p(i-1,j-1)) + g*H_g/(v_0*c_0))))/(2*(-v_0/(2*c_0))*(1/(a(i)^2))) ;
%         p_p(i,j)=(v_0/(2*c_0))*v_p(i,j)^(2)*(1/(a(i)^2))-g*H_g/(v_0*c_0);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    else
        v_p(i,j) = 0 ;
%         p_p(i,j) = v_p(i-1,j-1) + p_p(i-1,j-1) - landa*L/(2*D)*(v_0/c_0)*inc_t*v_p(i-1,j-1)*abs(v_p(i-1,j-1)) - v_p(i,j);
    end
        p_p(i,j) = v_p(i-1,j-1) + p_p(i-1,j-1) - landa*L/(2*D)*(v_0/c_0)*inc_t*v_p(i-1,j-1)*abs(v_p(i-1,j-1)) - v_p(i,j);
    for j = 2:lx-1
        b(1,1) = v_p(i-1,j-1) + p_p(i-1,j-1) - landa*L/(2*D)*(v_0/c_0)*inc_t*v_p(i-1,j-1)*abs(v_p(i-1,j-1));
        b(2,1) = v_p(i-1,j+1) - p_p(i-1,j+1) - landa*L/(2*D)*(v_0/c_0)*inc_t*v_p(i-1,j+1)*abs(v_p(i-1,j+1));
        x=A\b;
        v_p(i,j) = x(1);
        p_p(i,j) = x(2);
    end
end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tcoef = t_p(1:(t_c*N-1));

figure(5)
clf
plot(tcoef,ca)
xlabel('{\it t}','Fontsize',15)
ylabel('{\it ca}','Fontsize',15)

figure(6)
clf
plot(tcoef,cb)
xlabel('{\it t}','Fontsize',15)
ylabel('{\it cb}','Fontsize',15)


figure(7)
clf
plot(tcoef,cc)
xlabel('{\it t}','Fontsize',15)
ylabel('{\it cc}','Fontsize',15)
% 
% % Representación gráfica (animación)

%  m = size(v_p,1) ;
%  figure(1)
%  clf
% % hv = plot(x,v(1:n),'EraseMode','background');
%  hp = plot(x_p,p_p(1,:),'EraseMode','xor');
%  hold on
%  hv = plot(x_p,v_p(1,:),'r','EraseMode','xor');
%  axis([0 1 -1.2 1.2])
% % axis square
% % box on
%  set(gca,'xtick',[0 0.2 0.4 0.6 0.8 1])
%  set(gca,'ytick',[-1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1.0 1.2])
% % grid
%  xlabel('{\it x/L}','Fontsize',15)
%  ylabel('{\it v/v}_o ,   {\it p/p}_o','Fontsize',15)
%  title('Presión y velocidad a lo largo del conducto')
% 
%  for j=2:m
%    pause(0.01);
%    % set(hv,'XData',x,'YData',v( n*(j-1)+1 : n*j ))
%    set(hp,'YData',p_p(j,:))
%    set(hv,'YData',v_p(j,:))
%    drawnow
%  end

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rho_0 = 1000 ; % [kg/m3]
g = 9.81 ;     % [m/s2]
H_g = 400;     % [m]
p_a = 1*10^5 ; % [Pa]
c_0 = 1300 ; 

p_0 = rho_0*v_0*c_0;

alpha = asin(H_g/(Lgp+Ltf+Ltdif));
x = x_p*(Lgp+Ltf+Ltdif);
z = (Lgp+Ltf+Ltdif-x)*sin(alpha);
 
T = 20;                                            % [ºC] 
p_vapor = 1000*exp(16.573-3988.842/(273 + T - 39.47)) ; % [Pa]
% p_vapor_p = (p_vapor + rho_0*g*z - p_a - rho_0*g*H_g)/p_0;

for i = 1:lt
    for j = 1:lx
        p(i,j) = p_p(i,j)*p_0 - (rho_0*g*z(j) - p_a - rho_0*g*H_g) ;
        if (p(i,j)<p_vapor)
%         t = t_p(i)*(L/c_0);
        str = sprintf('\n Se ha alcanzado una presión inferior a la de vapor en el la posición i = %.0f, j = %.0f, correspondiente al instante t = %.0f s y a la posición x = %.0f m\n' , i, j, t_p(i)*(L/c_0),x_p(j)*L);
        disp(str)
        break
        end
    end
    if (p(i,j)<p_vapor)
        break
    end
end



% 
% figure(2)
% %  for i=1:lx
%     plot(p_vapor_p,x,p_x);
% %  end
%   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


landa = 0.02 ;     % coeficiente de fricción del conducto
D = 4.4 ;          % diámetro del conducto (m)
Agp = pi/4*D*D ;   % área del conducto a la entrada (m^2)
Atf = pi/4*D*D ;   % área del conducto a la entrada (m^2)
Atdif = 4.5*Atf ;    % área del conducto a la salida (m^2)
Ctdif = ( landa*Ltdif/D + 0.5 )/2/9.81/Atf/Atf ;
% v_E = v_p*v_0;
Q = v_p*v_0*Atf ;
%  H=-103,499 + 22,19·Q - 0,0785307·Q\u2\
B0= -103.499 ;
B1= 22.19 ;
B2= - 0.0785307 ;
HT = B0 + B1*Q + B2*Q.*Q ;

HS = Ctdif*Q.*Q;

HE = HS + HT ;
% z_E = 0 ;
z_E = H_g - (Lgp+Ltf)*sin(alpha);
p_E = (HE - z_E - ((v_p*v_0).*(v_p*v_0))/(2*g))*rho_0*g + p_a;
p_Ep = (p_E+rho_0*g*z_E-p_a-rho_0*g*H_g)/(p_0);

figure(2)
%  plot(t_p,p_p(:,12)','r')
%  hold on 
 plot(t_p,p_p(:,51)','r')
xlabel('{\it t/t}_i_v','Fontsize',15)
 ylabel('{\it p/p}_o','Fontsize',15) 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

