clear all
close all
clc
%% DATI GEOMETRIA MOTORE (motore quadro)


D = 68.8; %alesaggio cilindro [mm]
C = D;    % corsa del pistone

V = pi*(D^2/4)*C;    % cilindrata [cm^3]  
Vtot = V*4;   % cilindrata tot [cm^3] (4 cilindri)              

%% INTAKE VALVE GEOMETRY

dv = 0.40*D;                  %inner valve seat [mm]
dm = 1.10*dv;                 %head diameter [mm]
AB = dv +(sqrt(2)*dv);        %alesaggio cilindro tramite dv [mm]
ds = 0.20*dv;                 %stem diameter [mm]
beta_v = 45;                  %seat angle [deg]
r = 0.25*dv;                  %valve fitting radius [mm]
dh = 0.095*dv;                %da inner valve seat a valve head [mm]
w = ((dm-dv)/2)/cosd(beta_v); %seat width [mm]

%% GEOMETRIA CONDOTTO D'ASPIRAZIONE (entro la testata)

dc_in = 0.85*dv;  %diametro condotto inizio testata [mm]
dc_out = 0.93*dv; %diametro condotto prima della valvola [mm]
beta_c = 40;      %inclinazione condotto [deg]
L = 500;          %lunghezza condotto [mm]

%% CARATTERISTICHE DEL FLUSSO

a = 9.53*10^-1;                % coefficiente di Langen aria [kJ/kgK]
a1 = 6.65*10^-1;               % coefficiente di Langen aria [kJ/kgK]
b = 1.50*10^-4;                % coefficiente di Langen aria [kJ/kgK^2]
T_in = 298.15;                 % temperatura dell'aria in ingresso [K]
T_out = 293.08;                % temperatura dell'aria in uscita [k]
cp_in = (a+b*T_in)*1000;       % calore specifico a pressione costante [kJ/kgK]
cv_in = (a1+b*T_in)*1000;      % calore specifico a volume costante [kJ/kgK]
k = cp_in/cv_in;               % rapporto dei calori specifici
R = (cp_in-cv_in);             % costante dei gas [J/kgK]
cp_out=(a+b*T_out)*1000;
cv_out = (a1+b*T_out)*1000; 
k_out= cp_out/cv_out;                         % rapporto dei calori specifici
R_out = (cp_out-cv_out);                      % costante dei gas [J/kgK]
u0 = (2*(cp_in-cp_out)*(T_in-T_out))^0.5;     %velocità iniziale del flusso all'inlet
a_s = sqrt(k*R*T_in);                         % velocità caratteristica adiabatica [m/s]
n_giri = 6000;                                % regime di rotazione [rpm]
u_p = 2*C*(n_giri/60)*10^-3;                  %velocità media del pistone [m/s]
p0 = 95421.5;                                  % pressione in [Pa]
rho = p0/(R*T_in);                            % densità aria [kg/m^3]                      
Dp = 50*rho*(u_p^2)/2;                        % differenza di pressione tra ingresso e uscita [Pa]
p1 = p0-Dp;                               % pressione out [Pa]
beta = p0/p1;                                 % rapporto di compressione

#COEFFICIENTI PER IL CALCOLO DELLE PERDITE
                              
  frict_fact=0.015;
  lambda = 0.032; % conducibilità termica del condotto [W/m^2 K] 
  T_uscita_condotto = 350; %[K]
  
  
%% LAX WENDROFF  

%% geometria condtto espressa in metri (stessi diametri precedenti)

d_in = dc_in*10^-3; %diametro condotto_in [m]
d_out = dc_out*10^-3; %%diametro condotto_out [m]

%% B.C.
% pressione in ingresso e uscita definite precedentemente; velocità definita precedentemente
%% DISCRETIZZAZIONE

t = 1*10^-3;                %temporale [s]
s = L*10^-3;                %spaziale [m]
nx = 200;                   %nodi
Dx = s/nx;                  %passo
Dt = 0.05*(Dx/(a_s+u0));     %condizione Curant-Friederics-Lewis [s]
nt = ceil(t/Dt);            %numero di timestep
Dtx = Dt/Dx;
Courant = a_s*Dtx                %rapporto tempo/spazio di discretizzazione

%% INIZIALIZZAZIONE VETTORI

% Variabili            
  Rho = zeros(nx,nt);      
  u = zeros(nx,nt);                   
  P = zeros(nx,nt);                   
  T = zeros(nx,nt); 
  
  Ac = zeros(nx,1); 
  
% Vettore variabili conservate Q
  Q_1 = zeros(nx,nt); 
  Q_2 = zeros(nx,nt);
  Q_3 = zeros(nx,nt); 
  
% Vettore dei flussi F
  F_1 = zeros(nx,nt);
  F_2 = zeros(nx,nt);
  F_3 = zeros(nx,nt);
  
% Vettore termine sorgente S
 S_1 = zeros(nx,nt);
 S_2 = zeros(nx,nt);
 S_3 = zeros(nx,nt);
 
  
%% Variazione lineare sezione del codotto
 d=linspace(d_in,d_out,nx);
 Tw = linspace(T_out,T_uscita_condotto,nx); %Temperatura di parete del condotto [K]
for i=1:size(d,2)
    AA(i)=(pi*d(i)^2)/4;
end
Ac(nx,1)=AA(i);

%% START

for cicle = 1:1

for i=1:nt 

    %%inlet  
        Rho(1,i) = p0/(R*T_in);        
        u(1,i) = u0;        
        P(1,i) = p0;         
        T(1,i) = T_in;        
        Ac(1,1) = AA(1);        
        Q_1(1,i) = Rho(1,i);        
        Q_2(1,i) = Rho(1,i)*u(1,i);        
        Q_3(1,i) = 0.5*Rho(1,i)*u(1,i)^2 + P(1,i)/(k-1);      
        
        F_1(1,i) = Rho(1,i)*u(1,i);          
        F_2(1,i) = Rho(1,i)*u(1,i)^2 + P(1,i);         
        F_3(1,i) = u(1,i)*(0.5*Rho(1,i)*u(1,i)^2+P(1,i)*(k/(k-1))); 
        
        

     %%outlet
        Rho(nx,i) = p1/(R*T_uscita_condotto);         
        P(nx,i) = p1;        
        T(nx,i) = T_uscita_condotto;
        u(nx,i) = ((p0/(R*T_uscita_condotto))*AA(1)*u0)/(p1/(R*T_uscita_condotto)*AA(nx));        
        Ac(nx,1) = AA(nx);         
        Q_1(nx,i) = Rho(nx,i);         
        Q_2(nx,i) = Rho(nx,i)*u(nx,i);         
        Q_3(nx,i) = 0.5*Rho(nx,i)*u(nx,i)^2 + P(nx,i)/(k-1);         
        
        F_1(nx,i) = Rho(nx,i)*u(nx,i);         
        F_2(nx,i) = Rho(nx,i)*u(nx,i)^2 + P(nx,i);         
        F_3(nx,i) = u(nx,i)*(0.5*Rho(nx,i)*u(nx,i)^2+P(nx,i)*(k/(k-1))); 
        
        S_1(nx,i) = -Rho(nx,i)*u(nx,i)*(((Ac(nx,1)-Ac(nx-1,1))));
        S_2(nx,i) = -Rho(nx,i)*(u(nx,i)^2/Ac(nx,1))*(((Ac(nx,1)-Ac(nx-1,1))))-((frict_fact*pi*Rho(nx,i)*d(nx)*u(nx,i)**2)/(2*Ac(nx,1)));
        S_3(nx,i) = ((lambda*pi*d(nx)*(Tw(nx)-T(nx,i)))/((Ac(nx,1)-Ac(nx-1,1))))-u(nx,i)*(Rho(nx,i)*(u(nx,i)^2)/2+(k/(k-1))*P(nx,i))*(((Ac(nx,1)-Ac(nx-1,1))))/Ac(nx,1);
        
end

%% condotto
    for i=2:nx-1
        
        Rho(i,1) = p0/(R*Tw(i));         
        u(i,1) = u0;        
        P(i,1) = p0; 
        T(i,1) = Tw(i);        
        Ac(i,1)= AA(i); 
        
        Q_1(i,1) = Rho(i,1);         
        Q_2(i,1) = Rho(i,1)*u(i,1);         
        Q_3(i,1) = 0.5*Rho(i,1)*u(i,1)^2 + P(i,1)/(k-1);         
        
        F_1(i,1) = Rho(i,1)*u(i,1);        
        F_2(i,1) = Rho(i,1)*u(i,1)^2 + P(i,1);        
        F_3(i,1) = u(i,1)*(0.5*Rho(i,1)*u(i,1)^2+P(i,1)*(k/(k-1)));         
        
        S_1(i,1) = -Rho(i,1)*u(i,1)*(((Ac(i,1)-Ac(i-1,1))));
        S_2(i,1) = -Rho(i,1)*(u(i,1)^2/Ac(i,1))*(((Ac(i,1)-Ac(i-1,1))))-((frict_fact*pi*Rho(i,1)*d(i)*u(i,1)**2)/(2*Ac(i,1)));
        S_3(i,1) = ((lambda*pi*d(i)*(Tw(i)-T(i,1)))/((Ac(i,1)-Ac(i-1,1))))-u(i,1)*(Rho(i,1)*(u(i,1)^2)/2+(k/(k-1))*P(i,1))*(((Ac(i,1)-Ac(i-1,1))))/Ac(i,1);
           
    end
    
 %% cicle
    
    for n=2:nt 
        for i=2:nx-1
        
        %(i+1/2),(n+1/2)
        Q1_P_half = (1/2)*(Q_1(i,n-1) + Q_1(i+1,n-1)) - (Dt/(2*Dx))*(F_1(i+1,n-1)-F_1(i,n-1)) + (S_1(i+1,n-1)-S_1(i,n-1))*Dt/2;
        Q2_P_half = (1/2)*(Q_2(i,n-1) + Q_2(i+1,n-1)) - (Dt/(2*Dx))*(F_2(i+1,n-1)-F_2(i,n-1))+ (S_2(i+1,n-1)-S_2(i,n-1))*Dt/2;
        Q3_P_half = (1/2)*(Q_3(i,n-1) + Q_3(i+1,n-1)) - (Dt/(2*Dx))*(F_3(i+1,n-1)-F_3(i,n-1))+ (S_3(i+1,n-1)-S_3(i,n-1))*Dt/2;
        
        Rho_P_half = Q1_P_half;
        u_P_half = Q2_P_half/Rho_P_half;
        P_P_half = (Q3_P_half-0.5*Rho_P_half*u_P_half^2)*(k-1);
        
        F1_P_half = Q2_P_half;
        F2_P_half = ((Rho_P_half*u_P_half^2)/2)+P_P_half; 
        F3_P_half = u_P_half*(((Rho_P_half*u_P_half^2)/2)+(k*P_P_half/(k-1)));
        
        S_1_p1 = 0.5*(S_1(i+1,n-1)+S_1(i,n-1)); 
        S_2_p1 = 0.5*(S_2(i+1,n-1)+S_2(i,n-1)); 
        S_3_p1 = 0.5*(S_3(i+1,n-1)+S_3(i,n-1));
        
        
        %(i-1/2),(n+1/2)
        Q1_M_half = (1/2)*(Q_1(i,n-1) + Q_1(i-1,n-1)) - (Dt/(2*Dx))*(F_1(i,n-1)-F_1(i-1,n-1)) + (S_1(i,n-1)-S_1(i-1,n-1))*Dt/2;
        Q2_M_half = (1/2)*(Q_2(i,n-1) + Q_2(i-1,n-1)) - (Dt/(2*Dx))*(F_2(i,n-1)-F_2(i-1,n-1))+ (S_2(i,n-1)-S_2(i-1,n-1))*Dt/2;
        Q3_M_half = (1/2)*(Q_3(i,n-1) + Q_3(i-1,n-1)) - (Dt/(2*Dx))*(F_3(i,n-1)-F_3(i-1,n-1))+ (S_3(i,n-1)-S_3(i-1,n-1))*Dt/2;
        
        Rho_M_half = Q1_M_half;
        u_M_half= Q2_M_half/Rho_M_half;
        P_M_half = (Q3_M_half-0.5*Rho_M_half*u_M_half^2)*(k-1);
        
        F1_M_half = Q2_M_half;
        F2_M_half= ((Rho_M_half*u_M_half^2)/2)+P_M_half; 
        F3_M_half = u_M_half*(((Rho_M_half*u_M_half^2)/2)+(k*P_M_half/(k-1)));
        
        S_1_m1 = 0.5*(S_1(i,1)+S_1(i-1,1)); 
        S_2_m1 = 0.5*(S_2(i,1)+S_2(i-1,1));
        S_3_m1 = 0.5*(S_3(i,1)+S_3(i-1,1));
        
        
        % i
        Q_1(i,n) = Q_1(i,n-1) - Dtx*(F1_P_half - F1_M_half) +  (S_1_p1 - S_1_m1)*Dt;   
        Q_2(i,n) = Q_2(i,n-1) - Dtx*(F2_P_half - F2_M_half) +  (S_2_p1 - S_2_m1)*Dt;   
        Q_3(i,n) = Q_3(i,n-1) - Dtx*(F3_P_half - F3_M_half) +  (S_3_p1 - S_3_m1)*Dt;
        
##        Q_1(i,n) = 0.5*Courant*(1+Courant)*Q_1(i-1,n-1) + (1-Courant^2)*Q_1(i,n-1) - 0.5*Courant*(1-Courant)*Q_1(i+1,n-1)+  (S_1_p1 - S_1_m1)*Dt;   
##        Q_2(i,n) = 0.5*Courant*(1+Courant)*Q_2(i-1,n-1) + (1-Courant^2)*Q_2(i,n-1) - 0.5*Courant*(1-Courant)*Q_2(i+1,n-1)+  (S_2_p1 - S_2_m1)*Dt;   
##        Q_3(i,n) = 0.5*Courant*(1+Courant)*Q_3(i-1,n-1) + (1-Courant^2)*Q_3(i,n-1) - 0.5*Courant*(1-Courant)*Q_3(i+1,n-1)+  (S_3_p1 - S_3_m1)*Dt;
       
        Rho(i,n) = Q_1(i,n);
        u(i,n) = Q_2(i,n)/Rho(i,n);
        P(i,n) = (Q_3(i,n)-0.5*Rho(i,n)*u(i,n)^2)*(k-1);
    
        F_1(i,n) = Q_2(i,n);
        F_2(i,n) = Rho(i,n)*u(i,n)^2 + P(i,n);
        F_3(i,n) = u(i,n)*(k/(k-1)*P(i,n) + 0.5*Rho(i,n)*u(i,n)^2);
        
        S_1(i,n) = -Rho(i,n)*u(i,n)*(((Ac(i,1)-Ac(i-1,1))));
        S_2(i,n) = -Rho(i,n)*(u(i,n)^2/Ac(i,1))*(((Ac(i,1)-Ac(i-1,1))))-((frict_fact*pi*Rho(i,n)*d(i)*u(i,n)**2)/(2*Ac(i,1)));
        S_3(i,n) = ((lambda*pi*d(i)*(Tw(i)-T(i,n)))/((Ac(i,1)-Ac(i-1,1))))-u(i,n)*(Rho(i,n)*(u(i,n)^2)/2+(k/(k-1))*P(i,n))*(((Ac(i,1)-Ac(i-1,1))))/Ac(i,1); 
        
        end 
    end
end

%% P,u plots
figure(6)
plot(P(:,nt))
axis([0 nx+50 60000 200000])
xlabel('Nodi')
ylabel('Pressione [Pa]')
hold on

figure(7)
plot(u(:,nt))
axis([0 nx+50 1 20])
xlabel('Nodi')
ylabel('Velocità [m/s]')
hold on