Project 2 - Rankine cycle Simulator

%Rankine cycle

clear all
close all
clc

%Input parameters
%Considering the pump as the starting point, the inlet of Turbine will be
%point 3

% T3 = 673; %K; %400 C + 273 = 673 K, this is the Turbine inlet temperature
% P3 = 30; %bars; %Pressure entering the turbine
% P4 = 0.05; %bar; %Pressure exiting the turbine and entering the condenser. And the Pressure has reduced to a large extent which also indicates that the expansion occured in the turbine is very good/efficient.

fprintf("\n RANKINE CYCLE SIMULATOR");
fprintf("\n\n\n 1-2 represents Isentropic compression in the Pump");
fprintf("\n 2-3 represents Isobaric heat addition at the Boiler");
fprintf("\n 3-4 represents Isentropic expansion in the Turbine");
fprintf("\n 4-1 represents Isobaric heat removal at the condensor");

%To input the values
t3 = input("\n Enter Temperature at the Turbine Inlet (In Celcius)");
P3 = input("\n Enter the Pressure at the Turbine Inlet (In Bars)");
P4 = input("\n Enter the Pressure at the Turbine Outlet (In Bars)");

%Converting the temperature to Kelvin
%T3 = t3+273;

%Since we know the the Pressure at the inlet and outlet of the turbine, we
%can deduce P1 and P2 such that;
P1 = P4; %This is beacuse the heat is rejected at constant pressure in the condenser
P2 = P3; %Because the heat is added at constant pressure in the boiler

%To calculate the points
%At the inlet of the Pump
s1 = XSteam('sL_p', P1); %Saturated Liquid Entropy
h1 = XSteam('hL_p', P1); %Saturated Liquid Enthalpy
V1 = XSteam('vL_p', P1); %Saturated Liquid Volume
t1 = XSteam('Tsat_p', P1); %Saturated Temperature
%T1 = t1+273; %Converting into Kelvin

%Since Insentropic Compression occurs in the pump, we can say that s1 = s2;
%because Isentropic means the process is reversible and adiabatic. So what
%is the meaning of adiabatic: It is a process where neither heat is added or
%removed from the system.
s2 = s1;
%To calculate the enthaply at point 2 we need to know the work done. which
%can be described by the volume compressed by the pump, which is nothing
%but the work done by the pump.
WDp = V1*(P2-P1)*100;
h2 = WDp+h1;

%At the inlet of the boiler
%We need to define the Specific heat and Specific volume to know the state
%at point 2, this is used to define gamma, which is used to find the
%temperature at the point 2 using isentropic correlation.
Cp2 = XSteam('Cp_ps', P2, s2); %Specific isobaric heat capacity as a function of pressure and entropy.
Cv2 = XSteam('Cv_ps', P2, s2); %Specific isochoric heat capacity as a function of pressure and entropy.

gamma = Cp2/Cv2;
I = (gamma-1)/(gamma);

%The formula used for isentropic process is
%(T2/T1) = (P2/P1)^((gamma-1)/gamma);
t2 = t1*(P2/P1)^I;

%Now we can see at the inlet of the Turbine
%Tempearture and Pressure is known
h3 = XSteam('h_pT', P3, t3); %Entalpy as a function of pressure and temperature.
s3 = XSteam('s_pT', P3, t3); %Specific Entropy as a function of pressure and temperature

%Now at the inlet of the condenser
t4 = t1; %here heat removal occurs at const temp and const pr.
s4 = s3; %since it is isentropic expansion process

%To calculate the enthalpy change which results in the heat removal, we 1st
%need to the enthalpy and entropy of the state
sf4 = XSteam('sL_p', P4); %Saturated Liquid Entropy
sg4 = XSteam('sV_p', P4); %Saturated Vapour Entropy
hf4 = XSteam('hL_p', P4); %Saturated Liquid Enthalpy
hg4 = XSteam('hV_p', P4); %Saturated Vapour Entropy
%we need both liquid and vapour state here as during condensation there is
%a conversion of vapour to liquid

x4 = ((s4-sf4)/(sg4-sf4));

%This tells the heat removed from the system
h4 = hf4+(x4*(hg4-hf4));

%Now as all parameters are known we can get into the calculation of Work
%done, Heat supplied and removed, Efficiency and back work ratio and so on
WDt = h3-h4;

%Net work done;
W = WDt-WDp;

%Heat Supplied and Removed from the system
Qin = h3-h2;
Qout = h4-h1;

%This gives the thermal efficiency as we know that the total workdone by the heat input is the thermal efficiency
%It says how much heat has been utilized to get so much work as output
Ntherm = W/Qin;
Etherm = 1-(Qout/Qin); %This is also another way to find the thermal efficiency, as Work can be written as Qin-Qout

%Back work Ratio
%It is the ratio of the compressor work to the turbine work
B_WR = WDp/WDt;

%To find the Temp and entropy at point 5 and 6, that is at turbine inlet
%and outlet
t5 = XSteam('Tsat_p', P3);
%T5 = t5+273;

t6 = t5;

s5 = XSteam('sL_T', t5);
s6 = XSteam('sV_T', t6);

%Results
fprintf("\n Results:");
fprintf("\n State Point 1");
fprintf("\n P1 = %f", P1);
fprintf("\n T1 = %f", t1);
fprintf("\n s1 = %f", s1);
fprintf("\n h1 = %f", h1);

fprintf("\n\n\n State Point 2");
fprintf("\n P2 = %f", P2);
fprintf("\n T2 = %f", t2);
fprintf("\n s2 = %f", s2);
fprintf("\n h2 = %f", h2);

fprintf("\n\n\n State Point 3");
fprintf("\n P3 = %f", P3);
fprintf("\n T3 = %f", t3);
fprintf("\n s3 = %f", s3);
fprintf("\n h3 = %f", h3);

fprintf("\n\n\n State Point 4");
fprintf("\n P4 = %f", P4);
fprintf("\n T4 = %f", t4);
fprintf("\n s4 = %f", s4);
fprintf("\n h4 = %f", h4);

%Performance of the Rankine Cycle;
fprintf("\n\n\n Performance of the Rankine Cycle");
fprintf("\n Pump Work = %f", WDp);
fprintf("\n Turbine Work = %f", WDt);
fprintf("\n Net Work = %f", W);
fprintf("\n Back Work Ratio = %f", B_WR);
fprintf("\n Thermal Efficiency = %f", Ntherm);

%Plotting the graph

%T-S Diagram for Rankine Cycle
%For boiler heat is added at super heated condition
T = linspace(t6,t3,500); %Since we know that the line is not straight and is a curved line, therfore I use linspace
for i = 1:length(T)
s(i) = XSteam('s_pT', P3, T(i)); %Entropy as a function of pressure and temperature.
h(i) = XSteam('h_pT', P3, T(i)); %Entalpy as a function of pressure and temperature.
end

%Saturation Temperature Curve
t = linspace(1,t3,500);
for j=1:length(t)
sat_sL(j) = XSteam ('sL_T', t(j));
sat_sV(j) = XSteam ('sV_T', t(j));
sat_hL(j) = XSteam ('hL_T', t(j));
sat_hV(j) = XSteam ('hV_T', t(j));
end

%Plotting Saturation curves for liquid and vapour (The Liquid-Vapour dome)
plot(sat_sL,t,'color','r') %This is Saturation liquid temperature
hold on
plot(sat_sV,t,'color','r') %This is Saturation vapour temperature
%We are considering this because half the dome represents the liquid curve
%while the other half represnts the vapour state.

figure(1)
hold on
plot(s1, t1, 'o', 'color', 'g')
plot(s2, t2, 'o', 'color', 'g')
plot(s3, t3, 'o', 'color', 'g')
plot(s4, t4, 'o', 'color', 'g')
title('T-S Diagram of Rankine Cycle')
xlabel('Entropy (kJ/kg.K)')
ylabel('Temperature (K)')

line([s1 s2], [t1 t2], 'color', 'b', 'linewidth', 2)
line([s2 s5], [t2 t5], 'color', 'b', 'linewidth', 2)
line([s5 s6], [t5 t6], 'color', 'b', 'linewidth', 2)
plot(s,T, 'color', 'b', 'linewidth', 2)
line([s3 s4], [t3 t4], 'color', 'b', 'linewidth', 2)
line([s4 s1], [t4 t1], 'color', 'b', 'linewidth', 2)

figure
%plotting the HS curve
%Const pr heat rejection curve
S = linspace(s4, s1, 100);
for l = 1:length(S)
H(l) = XSteam('h_ps', P4, S(l));
end

Ent = linspace(s2, s3, 100);
for v = 1:length(Ent)
Enth(v) = XSteam('h_ps', P3, Ent(v));
end

%Plot for the saturation curves of liquid and vapour (The Liquid Vapour Dome)
plot(sat_sL, sat_hL,'color', 'r')
hold on
plot(sat_sV, sat_hV,'color', 'r')
%For HS liquid vapour dome, we need the saturation point of both the liquid
%and the vapour.

plot(s1, h1, 'o', 'color', 'g')
plot(s2, h2, 'o', 'color', 'g')
plot(s3, h3, 'o', 'color', 'g')
plot(s4, h4, 'o', 'color', 'g')
line([s1 s2], [h1 h2], 'color', 'b', 'linewidth', 2)
line([s3 s4], [h3 h4], 'color', 'b', 'linewidth', 2)
plot(S, H, 'color', 'b', 'linewidth', 2)
plot(Ent, Enth, 'color', 'b', 'linewidth', 2)

%Error obtained

%1. Since I took the temperature in K, the liquid vapour saturation curve was too small. Then I saw in the XSteam.m file that the temperature by default is taken in degree Celcius.

%2. I had this assumption because usually Turbine working points are calculated using Kelvin.

%Rakine cycle explaination

%1. I have assumed the cycle starting point to be at the pump, so the pump compresses [isentropic compression- which means reversible and adabatic process] (not much, as you can witness this from the result as well. Beacuse this is not a traditional compressor turbine), here temperature and Pr increases.

%2. The fluid then is let into the boiler where Isobaric (Const pr) heat addition occurs, which increases the temperature (a lot, because the external heat is added).

%3. Then there is an isentropic expansion occuring in the turbine. Which expands and produces energy. here as much as it can expand, it also reduces that much pressure and temperature.

%4. Then the fluid is let into the condersor where the Isobaric heat removal occurs. And the cycle continues. This cycles efficiency can further be increased using a reheat cycle. Which helps in reducing the heat supplied at the boiler. This is due to the two turbine (high temp and low temp Turbines). The heated fluid 1st enters the high pr turbine which expand extracting work, now the reduced temp and pr fluid enters the boiler again where there is an heat exchange occuring and then the fluid enters the low pr turbine, again extracting work.