Untitled

clc,clear,close all



global Q T_norm P_norm R cp cv Mw eta cost_energy hours Qm


Q = 500; %Nm^3/h
T_norm = 273.13; %  in K
P_norm = 101325; %  Pa
R = 8.314; %  J/(mol·K)
cp = 1006; %  J/(kg·K)
cv = 717.5; %  J/(kg·K)
Mw = 28.96; %  g/mol
eta = 0.72; % Efficienza
cost_energy = 0.15; % C €/kWh
hours = 6000; 


rho_norm = P_norm / (R * T_norm); 
Qm = Q * rho_norm / 1000 

lb = 1; %(in bar)
ub = 10; % (in bar)



P_sol = fmincon(@(Pint) obj_fn(Pint), 5, [], [], [], [], lb, ub, [], options)


Psuc1 = 1 * 10^5; 
Pdis1 = P_sol * 10^5; 
Tsuc1 = 298.15; 
[PwrComp1, Tdis1] = PowerC(Psuc1, Pdis1, Tsuc1);

Psuc2 = Pdis1; 
Pdis2 = 10 * 10^5; 
Tsuc2 = Tdis1;
[PwrComp2, Tdis2] = PowerC(Psuc2, Pdis2, Tsuc2);


Pwr_tot = PwrComp1 + PwrComp2
Savings = (Pwr_tot - 2 * PowerC(1 * 10^5, 10 * 10^5, 298.15)) * hours * cost_energy



function [PwrComp, Tdis] = PowerC(Psuc, Pdis, Tsuc)
    global cp cv Mw Qm eta
    
    
    Tdis = Tsuc * (Pdis / Psuc) ^ ((cp / cv) - 1) * (cp / cv);
    
   
    PwrComp = 2.31 * (cp / cv) * (Tdis - Tsuc) / Mw * Qm * eta / 1000;
end


function Pwr_tot = obj_fn(Pint)
    global cp cv Mw Qm
    
    
    Psuc1 = 1 * 10^5; 
    Pdis1 = Pint * 10^5; % in Pa
    Tsuc1 = 298.15; 
    
    [PwrComp1, Tdis1] = PowerC(Psuc1, Pdis1, Tsuc1);
    
    
    Psuc2 = Pdis1;
    Pdis2 = 10 * 10^5; 
    Tsuc2 = Tdis1; 
    
    [PwrComp2, Tdis2] = PowerC(Psuc2, Pdis2, Tsuc2);
    
    
    Pwr_tot = PwrComp1 + PwrComp2;
end