Project 2 - Rankine cycle Simulator

clear all
close all
clc

disp('RANKINE CYCLE SIMULATOR')
disp('1-2 is isentropic expansion')
disp('2-3 is constant pressure heat rejection')
disp('3-4 is isentropic compression')
disp('4-1 is constant pressure heat addition')

p1 = input('enter the pressure at the inlet of turbine =   ');
t1 = input('enter the temperature at the inlet of turbine =  ');
p2 = input('enter the pressure at the inlet of condenser =  ');

%State 1 conditions
h1 = XSteam('h_pT',p1,t1);
s1 = XSteam('s_pT',p1,t1);

%State 2 conditions
s2=s1;
sf2 = XSteam('sL_p',p2);
sg2 = XSteam('sV_p',p2);
x2 = (s2-sf2)/(sg2-sf2);
hf2 = XSteam('hL_p',p2);
hg2 = XSteam('hV_p',p2);
h2 = hf2 + (x2*(hg2-hf2));

%State 3 conditions 
p3=p2;
s3 = XSteam('sL_p',p3);
h3 = XSteam('hL_p',p3);

%State 4 conditions
s4=s3;
p4=p1;

%Work done by turbine
WT = h1-h2;

%Workdone to run pump
v3 = XSteam('vL_p',p3);
WP = v3*(p4-p3)*100;
h4 = h3+ WP;

%Net work done
WN = WT - WP;

%Heat into the system 
QIn = h1 - h4;

%Heat out of the system
QOut = h2 - h3;

%Thermal efficiency
Thermal_E = (WN/QIn)*100;

%Specific Steam Consumption
SSC = (3600/WN);

%Back work ratio
WB = (WP/WT);

%Temperature calculations at state points
t3 = XSteam('Tsat_p',p1);
t2=t3;
Cp4 = XSteam('Cp_ps',p4,s4);
Cv4 = XSteam('Cv_ps',p4,s4);
k = Cp4/Cv4;
t4 = t3/((p3/p4)^((k-1)/k));

%Temperature Entropy and Enthalpy calculations
t6 = XSteam('Tsat_p',p1);
s5 = XSteam('sL_p',p1);
s6 = XSteam('sV_p',p1);
h5 = XSteam('hL_p',p1);
h6  = XSteam('hV_p',p1);
t5 = t6;

%Plotting saturation curve
t = linspace(1,375,1000);
for i = 1:length(t)
    sf(i) = XSteam('sL_T',t(i));
    sg(i) = XSteam('sV_T',t(i));
    hf(i) = XSteam('hL_T',t(i));
    hg(i) = XSteam('hV_T',t(i));
end

%Plotting T-S curve
figure(1)
hold on 
plot(sf,t,'linewidth',2);
plot(sg,t,'linewidth',2);
plot([s1 s2],[t1 t2],'linewidth',2,'color','r');
plot([s2 s3],[t2 t3],'linewidth',2,'color','r');
plot([s3 s4],[t3 t4],'linewidth',2,'color','r');
plot([s4 s5],[t4 t5],'linewidth',2,'color','r');
plot([s5 s6],[t4 t5],'linewidth',2,'color','r');
sc1 = linspace(s6,s1,1000);
for k=1:length(sc1)
    tc1(1) = XSteam('T_ps',p1,sc1(1));
end
plot(sc1,tc1,'linewidth',2,'color','r')
title('TEMPERATURE VS ENTROPY');
xlabel('ENTROPY');
ylabel('TEMPERATURE');
legend('SATURATION CURVE')

%Plotting H-S curve
figure(2)
hold on 
plot(sf,hf,'linewidth',2);
plot(sg,hg,'linewidth',2);
plot([s1 s2],[h1 h2],'linewidth',2,'color','b');
plot([s2 s3],[h2 h3],'linewidth',2,'color','b');
plot([s3 s4],[h3 h4],'linewidth',2,'color','b');
hc = linspace(h4,h5,100);
for l=1:length(hc)
    sch(l) = XSteam('h_ps',p1,hc(l));
end
plot(sch,hc,'linewidth',2,'color','b')
plot([s5 s6],[h5 h6],'linewidth',2,'color','b');
for m=1:length(sc1)
    hc1(m) = XSteam('h_ps',p1,sc1(m));
end
plot(sc1,hc1,'linewidth',2,'color','b');
title('ENTHALPY VS ENTROPY');
xlabel('ENTROPY');
ylabel('ENTHALPY');
legend('SATURATION CURVE')

%Displaying outputs
%STATE ONE
fprintf('AT state point one:rn');
fprintf('P1 is : %f BARrn',p1);
fprintf('T1 is : %f Celciusrn',t1);
fprintf('H1 is : %f KJ/KGrn',h1);
fprintf('S1 is : %f KJ/KGrn',s1);
%STATE TWO
fprintf('AT state point two:rn');
fprintf('P2 is : %f BARrn',p2);
fprintf('T2 is : %f Celciusrn',t2);
fprintf('H2 is : %f KJ/KGrn',h2);
fprintf('S2 is : %f KJ/KGrn',s2);
fprintf('X is : %f rn',x2);
%STATE THREE
fprintf('AT state point one:rn');
fprintf('P3 is : %f BARrn',p3);
fprintf('T3 is : %f Celciusrn',t3);
fprintf('H3 is : %f KJ/KGrn',h3);
fprintf('S3 is : %f KJ/KGrn',s3);
%STATE FOUR
fprintf('AT state point four:rn');
fprintf('P4 is : %f BARrn',p4);
fprintf('T4 is : %f Celciusrn',t4);
fprintf('H4 is : %f KJ/KGrn',h4);
fprintf('S4 is : %f KJ/KGrn',s4);


%CALCULATIONS DISPLAY
fprintf('WP is %.3f KJ/KG',WP);
fprintf('WT is %.3f KJ/KG',WT);
fprintf('WN is %.3f KJ/KG',WN);
fprintf('THERMAL EFFICIENCY is %.3f PERCENT',Thermal_E);
fprintf('SSC is %.2f KG/KWH',SSC);
fprintf('BACK WORK RATIO IS %.3F',WB);