% sistema di condotti in serie
% Determinare la portata di acqua Q
% si trascurano le perdite accidentali
% Q=Q1=Q2=Q3 -> v2=(d1/d2)^2*v1   v3=(d1/d3)^2*v1
% delta_h=delta_h1+delta_h2+delta_h3
clear all;
close all;
clc;

g=9.81; % m/s^2 accelerazione di gravità

rho = 1000; % kg/m^3 densità acqua
nu = 1.02e-6; % m^2/s viscosità cinematica acqua

% pipe 1
L1=100; % m lunghezza
d1=8e-2; % m diametro
e1=0.24e-3; % m rugosità
ed1=e1/d1; % rugosità relativa

% pipe 2
L2=150; % m lunghezza
d2=6e-2; % m diametro
e2=0.12e-3; % m rugosità
ed2=e2/d2; % rugosità relativa

% pipe 3
L3=80; % m lunghezza
d3=4e-2; % m diametro
e3=0.20e-3; % m rugosità
ed3=e3/d3; % rugosità relativa

delta_p=150e3; % Pa salto di pressione tra ingresso e uscita
delta_z=5; % m dislivello tra ingresso e uscita

delta_h=delta_p/(rho*g)+delta_z;


% coeff d'attrito di primo tentativo (flusso in regime di rugosità piena)
f1 = friction_coeff_fully_rough_flow(ed1);
f2 = friction_coeff_fully_rough_flow(ed2);
f3 = friction_coeff_fully_rough_flow(ed3);

% delta_h=f1 * L1/d1 * u1^2/(2*g) + f2 * L2/d2 * u2^2/(2*g) +f3 * L3/d3 * u3^2/(2*g);
% delta_h=f1 * L1/d1 * u1^2/(2*g) + f2 * L2/d2 * (d1/d2)^4 * u1^2/(2*g) +f3 * L3/d3 * (d1/d3)^4 * u1^2/(2*g);
% delta_h=u1^2*(f1 * L1/d1 /(2*g) + f2 * L2/d2 * (d1/d2)^4 /(2*g) +f3 * L3/d3 * (d1/d3)^4 /(2*g));
% da cui si ricava l'espressione sottostante per il calcolo di u1:

u1 = sqrt(delta_h / (f1 * L1/d1 /(2*g) + f2 * L2/d2 * (d1/d2)^4 /(2*g)+f3 * L3/d3 * (d1/d3)^4 /(2*g)));

u2 = (d1/d2)^2 * u1;
u3 = (d1/d3)^2 * u1;

Re1 = u1*d1/nu;
Re2 = u2*d2/nu;
Re3 = u3*d3/nu;

Q = u1*pi*d1^2/4;

s=sprintf('f1=%0.6f f2=%0.6f f3=%0.6f   u1=%0.6f m/s u2=%0.6f m/s u3=%0.6f m/s  Re1=%0.0f Re2=%0.0f Re3=%0.0f   Q=%0.5f m^3/s',f1,f2,f3,u1,u2,u3,Re1,Re2,Re3,Q);
disp(s);

f1 = friction_coeff_Colebrook_formula(ed1,Re1);
f2 = friction_coeff_Colebrook_formula(ed2,Re2);
f3 = friction_coeff_Colebrook_formula(ed3,Re3);

u1_new = sqrt(delta_h / (f1 * L1/d1 /(2*g) + f2 * L2/d2 * (d1/d2)^4 /(2*g)+f3 * L3/d3 * (d1/d3)^4 /(2*g)));

eps=1e-3;
while abs(u1_new-u1)>eps
    u1=u1_new;
    
    u2 = (d1/d2)^2 * u1;
    u3 = (d1/d3)^2 * u1;
    
    Re1 = u1*d1/nu;
    Re2 = u2*d2/nu;
    Re3 = u3*d3/nu;
    
    Q = u1*pi*d1^2/4;
    
    s=sprintf('f1=%0.6f f2=%0.6f f3=%0.6f   u1=%0.6f m/s u2=%0.6f m/s u3=%0.6f m/s  Re1=%0.0f Re2=%0.0f Re3=%0.0f   Q=%0.5f m^3/s',f1,f2,f3,u1,u2,u3,Re1,Re2,Re3,Q);
    disp(s);

    f1 = friction_coeff_Colebrook_formula(ed1,Re1);
    f2 = friction_coeff_Colebrook_formula(ed2,Re2);
    f3 = friction_coeff_Colebrook_formula(ed3,Re3);
    
    u1_new = sqrt(delta_h / (f1 * L1/d1 /(2*g) + f2 * L2/d2 * (d1/d2)^4 /(2*g)+f3 * L3/d3 * (d1/d3)^4 /(2*g)));
end

Q = u1*pi*d1^2/4;
s=sprintf('\nQ=%0.3f m^3/s',Q);
disp(s);

% con la formula di Haaland
disp(' ');
disp('Con la formula di Haaland ');
% coeff d'attrito di primo tentativo (flusso in regime di rugosità piena)
f1 = friction_coeff_fully_rough_flow(ed1);
f2 = friction_coeff_fully_rough_flow(ed2);
f3 = friction_coeff_fully_rough_flow(ed3);

% delta_h=f1 * L1/d1 * u1^2/(2*g) + f2 * L2/d2 * u2^2/(2*g) +f3 * L3/d3 * u3^2/(2*g);
% delta_h=f1 * L1/d1 * u1^2/(2*g) + f2 * L2/d2 * (d1/d2)^4 * u1^2/(2*g) +f3 * L3/d3 * (d1/d3)^4 * u1^2/(2*g);
% delta_h=u1^2*(f1 * L1/d1 /(2*g) + f2 * L2/d2 * (d1/d2)^4 /(2*g) +f3 * L3/d3 * (d1/d3)^4 /(2*g));
% da cui si ricava l'espressione sottostante per il calcolo di u1:

u1 = sqrt(delta_h / (f1 * L1/d1 /(2*g) + f2 * L2/d2 * (d1/d2)^4 /(2*g)+f3 * L3/d3 * (d1/d3)^4 /(2*g)));

u2 = (d1/d2)^2 * u1;
u3 = (d1/d3)^2 * u1;

Re1 = u1*d1/nu;
Re2 = u2*d2/nu;
Re3 = u3*d3/nu;

Q = u1*pi*d1^2/4;

s=sprintf('f1=%0.6f f2=%0.6f f3=%0.6f   u1=%0.6f m/s u2=%0.6f m/s u3=%0.6f m/s  Re1=%0.0f Re2=%0.0f Re3=%0.0f   Q=%0.5f m^3/s',f1,f2,f3,u1,u2,u3,Re1,Re2,Re3,Q);
disp(s);

f1 = friction_coeff_Haaland_formula(ed1,Re1);
f2 = friction_coeff_Haaland_formula(ed2,Re2);
f3 = friction_coeff_Haaland_formula(ed3,Re3);

u1_new = sqrt(delta_h / (f1 * L1/d1 /(2*g) + f2 * L2/d2 * (d1/d2)^4 /(2*g)+f3 * L3/d3 * (d1/d3)^4 /(2*g)));

eps=1e-3;
while abs(u1_new-u1)>eps
    u1=u1_new;
    
    u2 = (d1/d2)^2 * u1;
    u3 = (d1/d3)^2 * u1;
    
    Re1 = u1*d1/nu;
    Re2 = u2*d2/nu;
    Re3 = u3*d3/nu;
    
    Q = u1*pi*d1^2/4;
    
    s=sprintf('f1=%0.6f f2=%0.6f f3=%0.6f   u1=%0.6f m/s u2=%0.6f m/s u3=%0.6f m/s  Re1=%0.0f Re2=%0.0f Re3=%0.0f   Q=%0.5f m^3/s',f1,f2,f3,u1,u2,u3,Re1,Re2,Re3,Q);
    disp(s);

    f1 = friction_coeff_Haaland_formula(ed1,Re1);
    f2 = friction_coeff_Haaland_formula(ed2,Re2);
    f3 = friction_coeff_Haaland_formula(ed3,Re3);
    
    u1_new = sqrt(delta_h / (f1 * L1/d1 /(2*g) + f2 * L2/d2 * (d1/d2)^4 /(2*g)+f3 * L3/d3 * (d1/d3)^4 /(2*g)));
end

Q = u1*pi*d1^2/4;
s=sprintf('\nQ=%0.3f m^3/s',Q);
disp(s);