% sistema di condotti in parallelo
% Determinare la portata di acqua Q
% si trascurano le perdite accidentali
% Q=Q1+Q2+Q3 
% 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)
% da cui si ricavano le espressioni per il calcolo di u1, u2, u3
u1=sqrt(delta_h / (f1*L1/d1/(2*g)));
u2=sqrt(delta_h / (f2*L2/d2/(2*g)));
u3=sqrt(delta_h / (f3*L3/d3/(2*g)));

Q = u1*pi*d1^2/4 + u2*pi*d2^2/4 + u3*pi*d3^2/4;

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

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)));
u2_new=sqrt(delta_h / (f2*L2/d2/(2*g)));
u3_new=sqrt(delta_h / (f3*L3/d3/(2*g)));

Q_new = u1_new*pi*d1^2/4 + u2_new*pi*d2^2/4 + u3_new*pi*d3^2/4;

eps=1e-4;
while abs(Q_new-Q)>eps
    Q=Q_new;
    u1=u1_new;
    u2=u2_new;
    u3=u3_new;
    
    Re1 = u1*d1/nu;
    Re2 = u2*d2/nu;
    Re3 = u3*d3/nu;
    
    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)));
    u2_new=sqrt(delta_h / (f2*L2/d2/(2*g)));
    u3_new=sqrt(delta_h / (f3*L3/d3/(2*g)));
    
    Q_new = u1_new*pi*d1^2/4 + u2_new*pi*d2^2/4 + u3_new*pi*d3^2/4;
    
end


s=sprintf('\nf1=%0.6f f2=%0.6f f3=%0.6f   u1=%0.6f m/s u2=%0.6f m/s u3=%0.6f m/s     Q=%0.5f m^3/s',f1,f2,f3,u1_new,u2_new,u3_new,Q_new);
disp(s);