% Función para el cálculo del caudal de fugas
function y=fuga1(p1,T1,V1,p2,m2,p3,m3,PE,TW,dif_ang,N,gamma,Mc,A,Ap,Vh,D,vp,Za,Cf)
Ro1=p1/((8314/Mc)*T1);
Ro2=m2/Vh;
Ro3=m3/Vh;
T2=p2/(Ro2*(8314/Mc));
T3=p3/(Ro3*(8314/Mc));
if p1>p2
if p2>p3
%Cálculo del caudal de masa
m2i=Cf*Ap*(p1*(2^0.5)/sqrt((8314/Mc)*T1))*sqrt((gamma/(gamma-1))*abs((p2/p1)^(2/gamma)-(p2/p1)^((gamma+1)/gamma)));
m2e=Cf*Ap*(p2*(2^0.5)/sqrt((8314/Mc)*T2))*sqrt((gamma/(gamma-1))*abs((p3/p2)^(2/gamma)-(p3/p2)^((gamma+1)/gamma)));
m3i=m2e;
m3e=Cf*Ap*(p3*(2^0.5)/sqrt((8314/Mc)*T3))*sqrt((gamma/(gamma-1))*abs((PE/p3)^(2/gamma)-(PE/p3)^((gamma+1)/gamma)));
else
m2i=Cf*Ap*(p1*(2^0.5)/sqrt((8314/Mc)*T1))*sqrt((gamma/(gamma-1))*abs((p2/p1)^(2/gamma)-(p2/p1)^((gamma+1)/gamma)))+Cf*Ap*(p3*(2^0.5)/sqrt((8314/Mc)*T3))*sqrt((gamma/(gamma-1))*abs((p2/p3)^(2/gamma)-(p2/p3)^((gamma+1)/gamma)));
m2e=0;
m3i=m2e;
m3e=Cf*Ap*(p3*(2^0.5)/sqrt((8314/Mc)*T3))*sqrt((gamma/(gamma-1))*abs((PE/p3)^(2/gamma)-(PE/p3)^((gamma+1)/gamma)));
end
else
if p2>p3
%Cálculo del caudal de masa
m2i=0;
m2e=Cf*Ap*(p2*(2^0.5)/sqrt((8314/Mc)*T2))*sqrt((gamma/(gamma-1))*abs((p3/p2)^(2/gamma)-(p3/p2)^((gamma+1)/gamma)))+Cf*Ap*(p2*(2^0.5)/sqrt((8314/Mc)*T2))*sqrt((gamma/(gamma-1))*abs((p1/p2)^(2/gamma)-(p1/p2)^((gamma+1)/gamma)));
m3i=Cf*Ap*(p2*(2^0.5)/sqrt((8314/Mc)*T2))*sqrt((gamma/(gamma-1))*abs((p3/p2)^(2/gamma)-(p3/p2)^((gamma+1)/gamma)));
m3e=Cf*Ap*(p3*(2^0.5)/sqrt((8314/Mc)*T3))*sqrt((gamma/(gamma-1))*abs((PE/p3)^(2/gamma)-(PE/p3)^((gamma+1)/gamma)));
else
m2i=Cf*Ap*(p3*(2^0.5)/sqrt((8314/Mc)*T3))*sqrt((gamma/(gamma-1))*abs((p2/p3)^(2/gamma)-(p2/p3)^((gamma+1)/gamma)));
m2e=Cf*Ap*(p2*(2^0.5)/sqrt((8314/Mc)*T2))*sqrt((gamma/(gamma-1))*abs((p1/p2)^(2/gamma)-(p1/p2)^((gamma+1)/gamma)));
m3i=0;
m3e=Cf*Ap*(p3*(2^0.5)/sqrt((8314/Mc)*T3))*sqrt((gamma/(gamma-1))*abs((PE/p3)^(2/gamma)-(PE/p3)^((gamma+1)/gamma)))+Cf*Ap*(p3*(2^0.5)/sqrt((8314/Mc)*T3))*sqrt((gamma/(gamma-1))*abs((p2/p3)^(2/gamma)-(p2/p3)^((gamma+1)/gamma)));
end
end
%Cálculo de las nuevas presiones
q2=heat_loses(T2,p2,TW,Vh,A,gamma,vp,D,0);
q3=heat_loses(T3,p3,TW,Vh,A,gamma,vp,D,0);
p22=p2+(gamma/Vh)*((m2i*dif_ang/N)*(p1/Ro1)-(m2e*dif_ang/N)*(p2/Ro2))-((gamma-1)/Vh)*1000*(dif_ang)*q2/N;
if Za<3
p33=p3;
else
p33=p3+(gamma/Vh)*((m3i*dif_ang/N)*(p2/Ro2)-(m3e*dif_ang/N)*(p3/Ro3))-((gamma-1)/Vh)*1000*(dif_ang)*q3/N;
end
%Cálculo de las masas
m2=m2+(m2i*dif_ang/N)-(m2e*dif_ang/N);
m3=m3+(m3i*dif_ang/N)-(m3e*dif_ang/N);
%Devolución de resultado
y=[m2i p22 p33 m2 m3 T2 T3];
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
