clear all
warning('off')
global R L1 L2 L3 D1 D2 D3 p1 p2 p3 T1 T2 T3 mu1 mu2 mu3
R = 286;
L1 = [0.6; 5];
L2 = [0.023; 0.02];
L3 = [0.02; 0.02];
D1 = [0.001; 0.001];
D2 = [0.0015; 0.0003];
D3 = [0.0003; 0.0003];
p1 = [1e5; 1e5];
p2 = [5.5*1e4; 5.5*1e4];
p3 = [4.8*1e4; 4.8*1e4];
T1 = [290; 290];
T2 = [420; 420];
T3 = [420; 290];
mu1 = [18*1e-6; 18*1e-6];
mu2 = [25*1e-6; 25*1e-6];
mu3 = [25*1e-6; 18*1e-6];
%Détermination de la pression de refoulement du tube de transfert
p2_adj = zeros(2,1);
options = optimset('TolX',1e-9,'MaxIter',1000000);
p2_adj(1) = fminsearch(@massflow_ffid,p2(1),options);
p2_adj(2) = fminsearch(@massflow_fid,p2(2),options);
%Détermination des débits
q1 = [0.001; 0.001];
q2 = [0.001; 0.001];
q3 = [0.001; 0.001];
for i = 1:10
    Re1 = 4*q1./(pi*mu1.*D1);
    Re2 = 4*q2./(pi*mu2.*D2);
    Re3 = 4*q3./(pi*mu3.*D3);
    f1 = friction_coef('lam',Re1);
    f2 = friction_coef('lam',Re2);
    f3 = friction_coef('lam',Re3);
    q1 = pi*D1.^2/4.*sqrt((p1.^2-p2_adj.^2)./(2*R*T1.*(2*f1.*L1./D1+log(p1./p2_adj))));
    q2 = pi*D2.^2/4.*sqrt((p2_adj.^2-p2.^2)./(2*R*T2(1).*(2*f2.*L2./D2+log(p2_adj./p2))));
    q3 = pi*D3.^2/4.*sqrt((p2.^2-p3.^2)./(2*R.*T3(1).*(2*f3.*L3./D3+log(p2./p3))));
end;
%Tracé de P, u, M en fonction de x/L
x_vect = [linspace(0,1,10);linspace(0,1,10)];
p1_vect = ((p1*ones(1,10)).^2 - (p1.^2-p2_adj.^2)*ones(1,10).*x_vect).^0.5;
p2_vect = ((p2_adj*ones(1,10)).^2 - (p2_adj.^2-p2.^2)*ones(1,10).*x_vect).^0.5;
p3_vect = ((p2*ones(1,10)).^2 - (p2.^2-p3.^2)*ones(1,10).*x_vect).^0.5;
u1_vect = q1.*R.*T1*ones(1,10)./(pi/4.*(D1*ones(1,10)).^2.*p1_vect);
u2_vect = q2.*R.*T2*ones(1,10)./(pi/4.*(D2*ones(1,10)).^2.*p2_vect);
u3_vect = q3.*R.*T3*ones(1,10)./(pi/4*(D3*ones(1,10)).^2.*p3_vect);
M1_vect = u1_vect./sqrt(1.4*R*T1*ones(1,10));
M2_vect = u2_vect./sqrt(1.4*R*T2*ones(1,10));
M3_vect = u3_vect./sqrt(1.4*R*T3*ones(1,10));
%DELAIS DE REPONSE
tau1 = 2.*D1./(f1.*u1_vect(:,1)).*(1 - p2_adj./p1 + R.*T1./(3*u1_vect(:,1).^2).*(1 - (p2_adj./p1).^3));
tau2 = 2.*D2./(f2.*u2_vect(:,1)).*(1 - p2./p2_adj + R.*T2./(3*u2_vect(:,1).^2).*(1 - (p2./p2_adj).^3));
tau3 = 2.*D3./(f3.*u3_vect(:,1)).*(1 - p3./p2 + R.*T2./(3*u2_vect(:,1).^2).*(1 - (p3./p2).^3));
%TEMPS DE REPONSE
tauc1_turb = 1.44*Re1.^(3/16).*L1./sqrt(R*T1).*sqrt((p1.^2 + p2_adj.^2)./(p1.^2 - p2_adj.^2));
tauc2_turb = 1.44*Re2.^(3/16).*L2./sqrt(R*T2).*sqrt((p2_adj.^2 + p2.^2)./(p2_adj.^2 - p2.^2));
tauc3_turb = 1.44*Re3.^(3/16).*L3./sqrt(R*T3).*sqrt((p2.^2 + p3.^2)./(p2.^2 - p3.^2));
tauc1_lam = 1.24*L1./sqrt(R*T1).*sqrt((p1.^2 + p2_adj.^2)./(p1.^2 - p2_adj.^2));
tauc2_lam = 1.24*L2./sqrt(R*T2).*sqrt((p2_adj.^2 + p2.^2)./(p2_adj.^2 - p2.^2));
tauc3_lam = 1.24*L3./sqrt(R*T3).*sqrt((p2.^2 + p3.^2)./(p2.^2 - p3.^2));
%AFFICHAGE
disp(['Délais Fast-FID : ', num2str(tau1(1)+tau2(1)+tau3(1)),' [s]']);
disp(['Délais FID : ', num2str(tau1(2)+tau2(2)),' [s]']);
disp(['Temps de réponse Fast-FID : ',num2str(tauc1_turb(1)+tauc2_turb(1)+tauc3_turb(1)),' [s]'])
disp(['Temps de réponse FID : ',num2str(tauc1_turb(2)+tauc2_turb(2)), ' [s]'])
figure(1)
subplot(3,1,1)
plot(x_vect,M1_vect(1,:),'b-',x_vect,M1_vect(2,:),'r-')
xlabel('Longueur normalisée (tube transfert) [-]')
ylabel('Nombre de Mach')
subplot(3,1,2)
plot(x_vect,M2_vect(1,:),'b-',x_vect,M2_vect(2,:),'r-')
xlabel('Longueur normalisée (tube tranquilisation) [-]')
ylabel('Nombre de Mach')
subplot(3,1,3)
plot(x_vect,M3_vect(1,:),'b-',x_vect,M3_vect(2,:),'r-')
xlabel('Longueur normalisée (tube connexion) [-]')
ylabel('Nombre de Mach')



