Project 2 - Rankine cycle Simulator

clear all
close all
clc

disp('Rankine cycle');
disp('1-2 is isentropic expansion in the turbine');
disp('2-3 is constant pressure heat removal in the condenser');
disp('3-4 is isentropic compression in the pump');
disp('4-1 is constant pressure heat addition in the boiler')
fprintf('\n');
%Obtaining the input data
p1=input('Pressure at turbine inlet,p1(in bar)=');
T1=input('Temperature at turbine inlet,T1(in degree.C)=');
p2=input('Temperature at condenser inlet,p2(in bar)=');

%Extraxting data and defining variables
%State 1
h1=XSteam('h_pT',p1,T1);
s1=XSteam('s_pT',p1,T1);

%State 2
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
p3=p2;
s3=XSteam('sL_p',p3);
h3=XSteam('hL_p',p3);

%State 4
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 
wnet=wt-wp;

%Heat into the system
Qin=h1-h4;

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

%Thermal efficiency
eff=(wnet/Qin)*100;

%Specific steam consumption
%SSC=mass flow rate/power
SSC=(3600/wnet);

%Calculating temperature at the state points
T3=XSteam('Tsat_p',p3);
T2=T3;
Cp4=XSteam('Cp_ps',p4,s4);
Cv4=XSteam('Cv_ps',p4,s4);
k=Cp4/Cv4;
T4=T3/(p3/p4)^((k-1)/k);

%Calculating temperature,entropy and enthalpy at the vapour and liquid
%line for plotting
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 the 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 diagram
figure(1)
hold on
plot(sf,T,'r','linewidth',1)
plot(sg,T,'r','linewidth',1)
plot([s1 s2],[T1 T2],'k','linewidth',2)
plot([s2 s3],[T2 T3],'k','linewidth',2)
plot([s3 s4],[T3 T4],'k','linewidth',2)
plot([s4 s5],[T4 T5],'k','linewidth',2)
plot([s5 s6],[T5 T6],'k','linewidth',2)
plot([s6 s1],[T6 T1],'k','linewidth',2)
sc1=linspace(s6,s1,1000);
for i=1:length(sc1)
    Tc1(1)=XSteam('T_ps',p1,sc1(1));
end
plot(sc1,Tc1,'k','linewidth',2)
text(s1+0.1,T1,'1')
text(s2+0.1,T2,'2')
text(s3,T3-20,'3')
text(s4,T4+20,'4')
title('Temp Vs Entropy plot')
xlabel('Entropy,S(Kj/kgk')
ylabel('Temperature (degree c)')
legend('Saturation curve','location','northwest')

%Plotting h-s diagram
figure(2)
hold on
plot(sf,hf,'r','linewidth',2)
plot(sg,hg,'r','linewidth',2)
plot([s1 s2],[h1 h2],'k','linewidth',2)
plot([s2 s3],[h2 h3],'k','linewidth',2)
plot([s3 s4],[h3 h4],'k','linewidth',2)
hc=linspace(h4,h5,100);
for k=1:length(hc)
    sch(k)=XSteam('s_ph',p1,hc(k));
end
plot(sch,hc,'k','linewidth',2)
plot([s5 s6],[h5 h6],'k','linewidth',2)
for m=1:length(sc1)
    hc1(m)=XSteam('h_ps',p1,sc1(m));
end
plot(sc1,hc1,'k','linewidth',2)
text(s1+0.1,h1,'1')
text(s2+0.1,h2,'2')
text(s3,h3-100,'3')
text(s4,h4+200,'4')
title('Enthalpy Vs Entropy plot')
xlabel('Entropy, S(KJ/kgk)')
ylabel('Enthalpy,H(KJ/kg)')
legend('Saturation curve','location','northwest')

%Displaying Results
fprintf('\n');
disp('RESULTS');
fprintf('At state point 1\n');
fprintf('P1:%8.3f bar\n',p1);
fprintf('T1:%8.3f C\n',T1);
fprintf('h1:%8.3f KJ/kg\n',h1);
fprintf('s1:%8.3f KJ/kgKn\n',s1);
fprintf('\n');
fprintf('At State point 2\n');
fprintf('p2:%8.3f bar\n',p2);
fprintf('T2:%8.3f C\n',T2);
fprintf('h2:%8.3f KJ/kg\n',h2);
fprintf('s2:%8.3f KJ/KgKn\n',s2);
fprintf('x2:%8.3fn\n',x2);
fprintf('\n');
fprintf('At state point 3\n');
fprintf('p3:%8.3f bar\n',p3);
fprintf('T3:%8.3f K\n',T3);
fprintf('h3:%8.3f KJ/kg\n',h3);
fprintf('s3:%8.3f KJ/kgKn\n',s3);
fprintf('\n');
fprintf('At State point 4\n');
fprintf('p4:%8.3f bar\n',p4);
fprintf('T4:%8.3f K\n',T4);
fprintf('h4:%8.3f K\n',h4);
fprintf('s4:%8.3f KJ/kgKn\n',s4);
fprintf('\n');
fprintf('Turbine work:%8.3f K\n',wt);
fprintf('Pump work:%8.3f KJ/kg\n',wp);
fprintf('Wnet:%8.3f KJ/kg\n',wnet);
fprintf('Efficiency:%8.3f Percent \n',eff);
fprintf('S.S.C:%8.3f KJ/kWh',SSC);