clear all

% This is a single-screw expander model, built using the relations from
% work by A. Giuffrida "Improving the semi-empirical modelling of a single-screw expander for
% small organic Rankine cycles" (2017). To run the model, it is neccesary
% to install to CoolProp and Refprop software and link them to MATLAB via Python. 
% Once installed, simply specify the T_e, T_a and BVR and
% select the range of V_sw and pressure ratios. The results of
% interest will then be displayed in table T. The guessed values of the flow rates and shell
% temperatures may need to be calibrated first. 
% IF You wish to obtai nresults for Variant 2, please remember to set 
% n = 1500 rpm and a_leak1 = 0


%%% Call Refprop %%%
global RP
global MOLSI
global iMass
global iFlag

RP = py.ctREFPROP.ctREFPROP.REFPROPFunctionLibrary('C:\Program Files (x86)\REFPROP');
MOLSI = RP.GETENUMdll(int8(0),'SI WITH C').iEnum;
iMass = int8(1); % 0: molar fractions; 1: mass fractions
iFlag = int8(0); % 0: don't call SATSPLN; 1: call SATSPL
z = {1.0}; % mole fractions, here a pure fluid, so mole fraction of 1.0

%%% Pressures and admission and ideal emission points for the base cycle %%%
T_e = 200; % C
T_a = 450; % C
Power = 10; % kW
eta_pump = 0.48;
p_e = py.CoolProp.CoolProp.PropsSI('P','T',T_e+273.15,'Q',1,'Water')/1000000; % MPa
epsilon = 1.5;
epsilon_max = 1.5;
epsilon_count = (epsilon_max - epsilon)*2 + 1;
BVR = 1.5;
% Input parameters %
        

A_leak_nom = 17; %mm2
a_leak1 = 7.6;%mm2/MPa
A_in = 93; %mm2
K = 1.5;
f_loss0 = 103.2 * 10^(-6);
f_loss1 = 3.03 * 10^(-4);
b_hl = 1.32; %W*K^(-5/4)
b_ra = 3.14 * 10^(-8);

t_amb = 20;

kappa = 1.33;
n = 3000; %RPM
eta_gen = 0.912;

tab_index = 0;

for j=1:epsilon_count
    
    V_sw_cm = 1100; %cm3
    V_sw_max_cm = 1100; %cm3
  
    count = (V_sw_max_cm - V_sw_cm)/100 + 1;
    

    for i=1:count
        t_shell = T_a - 40 ; %C guess
        m_dot_nom = 0.01; %kg/s (guess)
        V_sw = V_sw_cm/1000000; %m3
        p_a = p_e*epsilon;
        
        % Pressures
        p1 = p_e;
        p2s = p_a;
        p2 = p_a;
        p_sl = p_a;
        p_sv = p_a;
        p3 = p_a;
        p4s = p_e;
        p4 = p_e;
        
        % Point 3
        t3 = T_a;
        h3 = py.CoolProp.CoolProp.PropsSI('H','T',t3+273.15,'P',p3*1000000,'Water')/1000; % kJ/kg
        s3 = py.CoolProp.CoolProp.PropsSI('S','T',t3+273.15,'P',p3*1000000,'Water')/1000; % kJ/kgK
        cp3 = py.CoolProp.CoolProp.PropsSI('C','T',t3+273.15,'P',p3*1000000,'Water'); % J/kgK
        
        % Point 4s
        s4s = s3;
        h4s = py.CoolProp.CoolProp.PropsSI('H','S',s3*1000,'P',p4*1000000,'Water')/1000; % kJ/kg
        x4s = py.CoolProp.CoolProp.PropsSI('Q','S',s3*1000,'P',p4*1000000,'Water'); 
        
        
        %%% Expander model %%%
        
        % 0 -> su is adiabatic throttling at the expander inlet
        % su -> su1 is isobaric cooling of the WF by the imaginary shell
        % su1 -> in models the leakage loss
        % in -> int is the isentropic expansion
        % int -> ex2 is the isochoric expansion
        % ex2 -> ex1 is the inflow of the leaked fluid
        % ex1 -> ex is the isobaric heating of the WF by the imaginary shell
        
        error_shell = 10;
        
        m_leak = 0.001; %kg/s (guess)
        m_dot = m_dot_nom + m_leak;
        
        error_nom = 1;
        error_leak = 1;
        
        while abs(error_nom) > 0.01 || abs(error_leak) > 0.01 || abs(error_shell) > 0.025
        
            % Point 0
            t_0 = T_a;
            p_0 = p_a;
            h_0 = py.CoolProp.CoolProp.PropsSI('H','T',t_0 + 273.15,'P',p_0*1000000,'Water')/1000;
            s_0 = py.CoolProp.CoolProp.PropsSI('S','T',t_0 + 273.15,'P',p_0*1000000,'Water')/1000;
            m_0 = m_dot;
        
            % Point su
            p_su = p_0 - 0.1*p_a; % guess
            s_su_s = s_0;
            error_supply = 1;
        
             while abs(error_supply) > 0.01
              
                d_su_s = py.CoolProp.CoolProp.PropsSI('D','S',s_su_s*1000,'P',p_su*1000000,'Water');
                h_su_s = py.CoolProp.CoolProp.PropsSI('H','S',s_su_s*1000,'P',p_su*1000000,'Water')/1000;
                m_dot_calc = d_su_s*(A_in*10^(-6))*sqrt(2*1000*(h_0 - h_su_s));
                error_supply = m_dot - m_dot_calc;
        
                if error_supply < 0
                    p_su = p_su + 0.001;
                elseif error_supply > 0
                    p_su = p_su - 0.001;
                else
                    p_su = p_su;
                end

            end
        
            
            h_su = h_0;
            t_su = py.CoolProp.CoolProp.PropsSI('T','H',h_su*1000,'P',p_su*1000000,'Water') - 273.15;
            s_su = py.CoolProp.CoolProp.PropsSI('S','H',h_su*1000,'P',p_su*1000000,'Water')/1000;
            m_su = m_dot;
            
            % Point su1
        
            lambda_su = py.CoolProp.CoolProp.PropsSI('L','H',h_su*1000,'P',p_su*1000000,'Water');
            vis_su = py.CoolProp.CoolProp.PropsSI('V','H',h_su*1000,'P',p_su*1000000,'Water');
            cp_su = py.CoolProp.CoolProp.PropsSI('C','H',h_su*1000,'P',p_su*1000000,'Water');
            prandtl_su = cp_su*vis_su/lambda_su;
        
            AU_in = K*lambda_su*((m_dot/vis_su)^0.8)*prandtl_su^0.3;
        
            Q_su = (1-exp((-AU_in)/(m_dot*cp_su)))*m_dot*cp_su*(t_su-t_shell)/1000; % kW
            q_su = Q_su/m_dot; % kJ/kg
            
            h_su1 = h_su - q_su;
            p_su1 = p_su;
            s_su1 = py.CoolProp.CoolProp.PropsSI('S','H',h_su1*1000,'P',p_su1*1000000,'Water')/1000; % kJ/kgK
            t_su1 = py.CoolProp.CoolProp.PropsSI('T','H',h_su1*1000,'P',p_su1*1000000,'Water') - 273.15;
            d_su1 = py.CoolProp.CoolProp.PropsSI('D','H',h_su1*1000,'P',p_su1*1000000,'Water');
            v_su1 = 1/d_su1;
            m_su1 = m_dot;
            
            % Point in
            
            h_in = h_su1;
            t_in = t_su1;
            p_in = p_su1;
            s_in = s_su1;
            d_in = d_su1;
            v_in = v_su1;
            m_in = m_dot_nom;
            
            % Point int
            s_int = s_in;
            v_int = v_in*BVR;
            d_int = 1/v_int;
        
            r = RP.REFPROPdll('Water','SD','P',MOLSI,iMass,iFlag,s_int,d_int,z);
            o = double(r.Output);
            p_int = o(1);
        
            t_int = py.CoolProp.CoolProp.PropsSI('T','S',s_int*1000,'P',p_int*1000000,'Water') - 273.15;
            h_int = py.CoolProp.CoolProp.PropsSI('H','S',s_int*1000,'P',p_int*1000000,'Water')/1000;
            m_int = m_dot_nom;
            
            % Point ex2
            p_ex2 = p_e;
            h_ex2 = h_int-v_int*1000*(p_int-p_ex2);
            d_ex2 = py.CoolProp.CoolProp.PropsSI('D','P',p_ex2*1000000,'H',h_ex2*1000,'Water');
            t_ex2 = py.CoolProp.CoolProp.PropsSI('T','P',p_ex2*1000000,'H',h_ex2*1000,'Water') - 273.15;
            s_ex2 = py.CoolProp.CoolProp.PropsSI('S','P',p_ex2*1000000,'H',h_ex2*1000,'Water')/1000;
            v_ex2 = 1/d_ex2;
            m_ex2 = m_dot_nom;
            x_ex_2 = py.CoolProp.CoolProp.PropsSI('Q','P',p_ex2*1000000,'H',h_ex2*1000,'Water');
            
            % m_dot and m_leak calculation and iteration
            mep = d_int*(h_in - h_ex2)/1000;
            p_load = p_e + mep;
            m_dot_nom_new = d_in*V_sw*n/(BVR*60);
            p_crit = p_in*(2/(kappa+1))^(kappa/(kappa-1));
            
            if p_crit > p_ex2
                p_leak = p_crit;
            else
                p_leak = p_ex2;
            end
        
            s_leak = s_su1;
            d_leak = py.CoolProp.CoolProp.PropsSI('D','P',p_leak*1000000,'S',s_su1*1000,'Water');
            h_leak = py.CoolProp.CoolProp.PropsSI('H','P',p_leak*1000000,'S',s_su1*1000,'Water')/1000;
            
            A_leak = A_leak_nom + a_leak1*p_load;
        
            m_leak_new = d_leak*A_leak/1000000*sqrt(2000*(h_su1-h_leak));
            m_dot_new = m_dot_nom_new + m_leak_new;
        
            error_nom = m_dot_nom - m_dot_nom_new;
            error_leak = m_leak - m_leak_new;
            m_dot_nom = m_dot_nom_new;
            m_leak = m_leak_new;
            m_dot = m_dot_new;
            
            
            % Point ex1
            p_ex1 = p_ex2;
            h_ex1 = (m_dot_nom*h_ex2 + m_leak*h_su1)/m_dot;
            t_ex1 = py.CoolProp.CoolProp.PropsSI('T','P',p_ex1*1000000,'H',h_ex1*1000,'Water') - 273.15;
            s_ex1 = py.CoolProp.CoolProp.PropsSI('S','P',p_ex1*1000000,'H',h_ex1*1000,'Water')/1000;
            d_ex1 = py.CoolProp.CoolProp.PropsSI('D','P',p_ex1*1000000,'H',h_ex1*1000,'Water');
            m_ex1 = m_dot; 
            cp_ex1 = py.CoolProp.CoolProp.PropsSI('C','P',p_ex1*1000000,'H',h_ex1*1000,'Water'); %J/kgK
            
            % Point ex
            lambda_ex1 = py.CoolProp.CoolProp.PropsSI('L','H',h_ex1*1000,'P',p_ex1*1000000,'Water');
            vis_ex1 = py.CoolProp.CoolProp.PropsSI('V','H',h_ex1*1000,'P',p_ex1*1000000,'Water');
            prandtl_ex1 = cp_ex1*vis_ex1/lambda_ex1;
        
            AU_ex = K*lambda_ex1*((m_dot/vis_ex1)^0.8)*prandtl_ex1^0.4;
        
            Q_ex = (1-exp((-AU_ex)/(m_dot*cp_ex1)))*m_dot*cp_ex1*(t_shell-t_ex1)/1000; % kW
            q_ex = Q_ex / m_dot; %kJ/kg
            
            h_ex = h_ex1 + q_ex;
            p_ex = p_ex2;
            s_ex = py.CoolProp.CoolProp.PropsSI('S','H',h_ex*1000,'P',p_ex*1000000,'Water')/1000; % kJ/kgK
            t_ex = py.CoolProp.CoolProp.PropsSI('T','H',h_ex*1000,'P',p_ex*1000000,'Water')-273.15; % C
            x_ex = py.CoolProp.CoolProp.PropsSI('Q','H',h_ex*1000,'P',p_ex*1000000,'Water');
            
            % Internal and shaft power, losses, shell temperature iteration 
            P_in = d_int*(h_in - h_ex2)*V_sw*n/60; % %kW
            P_loss = (f_loss0 + f_loss1*n/(60*p_load*1000))*p_load*1000*n/60; %kW
            P_sh = P_in - P_loss;
            Q_amb_balance = P_loss + (1 - eta_gen)*P_sh + Q_su - Q_ex;
            Q_amb_calc = (b_hl*(t_shell - t_amb)^(1.25) + b_ra*((t_shell + 273.15)^4 - (t_amb + 273.15)^4))/1000;
            error_shell = Q_amb_balance - Q_amb_calc;
        
            if error_shell > 0
        
                if error_shell > 10
                    t_shell = t_shell + 10;
                elseif error_shell > 1
                    t_shell = t_shell + 1;
                else 
                    t_shell = t_shell + 0.1;
                end
        
            else
        
                if t_shell < (-1)
                    t_shell = t_shell - 1;
                else
                    t_shell = t_shell - 0.1;
                end
        
            end
        end
        
        %%% Shaft power and expander efficiency %%%
        P_el_id = m_dot*(h3-h4s);
        eta_exp = P_sh/P_el_id;
        
        %%% Rest of the base cycle - points 1, 2, 2s, SV, SL and 4 %%%
        % Point 4
        h4 = h3-(h3-h4s)*eta_exp;
        s4 = py.CoolProp.CoolProp.PropsSI('S','H',h4*1000,'P',p4*1000000,'Water')/1000; % kJ/kgK
        x4 = py.CoolProp.CoolProp.PropsSI('Q','S',s4*1000,'P',p4*1000000,'Water'); 
        t4 = py.CoolProp.CoolProp.PropsSI('T','S',s4*1000,'P',p4*1000000,'Water') - 273.15; 
        
        % Point 1
        x1 = 0;
        t1 = py.CoolProp.CoolProp.PropsSI('T','Q',x1,'P',p1*1000000,'Water') - 273.15; 
        h1 = py.CoolProp.CoolProp.PropsSI('H','Q',x1,'P',p1*1000000,'Water')/1000; 
        s1 = py.CoolProp.CoolProp.PropsSI('S','Q',x1,'P',p1*1000000,'Water')/1000; 
        
        % Point 2s
        s2s = s1;
        h2s = py.CoolProp.CoolProp.PropsSI('H','S',s2s*1000,'P',p2*1000000,'Water')/1000; % kJ/kg
        t2s = py.CoolProp.CoolProp.PropsSI('T','S',s2s*1000,'P',p2*1000000,'Water') -273.15; 
        
        % Point 2
        h2 = h1+(h2s-h1)/eta_pump;
        s2 = py.CoolProp.CoolProp.PropsSI('S','H',h2*1000,'P',p2*1000000,'Water')/1000; % kJ/kgK
        t2 = py.CoolProp.CoolProp.PropsSI('T','S',s2*1000,'P',p2*1000000,'Water')-273.15; 
        
        % Point SL
        x_sl = 0;
        t_sl = py.CoolProp.CoolProp.PropsSI('T','Q',x_sl,'P',p_sl*1000000,'Water') - 273.15; 
        h_sl = py.CoolProp.CoolProp.PropsSI('H','Q',x_sl,'P',p_sl*1000000,'Water')/1000; 
        s_sl = py.CoolProp.CoolProp.PropsSI('S','Q',x_sl,'P',p_sl*1000000,'Water')/1000; 
        
        % Point SV
        x_sv = 1;
        t_sv = py.CoolProp.CoolProp.PropsSI('T','Q',x_sv,'P',p_sv*1000000,'Water')-273.15; 
        h_sv = py.CoolProp.CoolProp.PropsSI('H','Q',x_sv,'P',p_sv*1000000,'Water')/1000; 
        s_sv = py.CoolProp.CoolProp.PropsSI('S','Q',x_sv,'P',p_sv*1000000,'Water')/1000; 
        
        Q_out = m_dot*(h_ex-h1);
        Q_in = m_dot*(h3-h2);
        W_k = m_dot*(h2 - h1);
        P_e = m_dot*(h3-h4); %control check 
        P_el = P_e*eta_gen;
        P_c = P_el - W_k;
        eta_el = P_c/Q_in;
        
        tab_index = tab_index + 1;
        A(tab_index) = V_sw_cm;
        B(tab_index) = eta_exp*100;
        C(tab_index) = P_c;
        D(tab_index) = eta_el*100;
        E(tab_index) = Q_out;
        F(tab_index) = x4;
        G(tab_index) = t_ex;
        H(tab_index) = p_load;
        I(tab_index) = m_leak;
        J(tab_index) = m_dot;
        

        V_sw_cm = V_sw_cm + 100;
        X(tab_index) = epsilon;
    
    end
    epsilon = epsilon + 0.5;
end

T = table(X.',A.',B.',C.',D.',E.',F.',G.',H.',I.',J.','VariableNames',{'epsilon','V_sw','Eta_exp','Power','eta_el','HeatOut','Quality','ExhaustTemp','p_load','m_leak','m_dot'});
 disp(T)

