Project 2 - Rankine cycle Simulator

Aim: Create a Rankine Cycle Simulator using MATLAB.

Objective:

• Calculate each state points of the Rankine cycle (any type of your choice) based on user inputs.
• Plot the corresponding T-S and H-S plots for the given set of inputs.
• Calculate Efficiency.

Theory: The Rankine cycle is an idealized thermodynamic cycle of a heat engine that converts heat into mechanical work while undergoing a phase change. It is an idealized cycle in which friction losses in each of the four components are neglected. The heat is supplied externally to a closed loop, which usually uses water as the working fluid.

Four components of the Rankine cycle:

1. Pump
2. Boiler
3. Turbine
4. Condenser.

Four processes of the Rankine cycle:

• Isentropic compression: (compression in centrifugal pumps) – The liquid condensate is compressed adiabatically from state 1 to state 2 by centrifugal pumps (usually by condensate pumps and then by feedwater pumps). The liquid condensates water pumped from the condenser into the higher pressure boiler. In this process, the surroundings do work on the fluid, increasing its enthalpy (h = u+pv) and compressing it (increasing its pressure). On the other hand the entropy remains unchanged. The work required for the compressor is given by W_Pumps = h₂ – h₁.

• Isobaric heat addition: (in a heat exchanger – boiler) – In this phase (between state 2 and state 3) there is a constant pressure heat transfer to the liquid condensate from an external source, since the chamber is open to flow in and out. The feedwater (secondary circuit) is heated from to the boiling point (2 → 3a) of that fluid and then evaporated in the boiler (3a → 3). The net heat added is given by Q_add = h₃ – h_2.

• Isentropic expansion: (expansion in a steam turbine) – Steam from the boiler expands adiabatically from state 3 to state 4 in a steam turbine to produce work and then is discharged to the condenser (partially condensed). The steam does work on the surroundings (blades of the turbine) and loses an amount of enthalpy equal to the work that leaves the system. The work done by turbine is given by W_T = h₄ – h₃. Again the entropy remains unchanged.

• Isobaric heat rejection: (in a heat exchanger) - In this phase the cycle completes by a constant-pressure process in which heat is rejected from the partially condensed steam. There is heat transfer from the vapor to cooling water flowing in a cooling circuit. The vapor condenses and the temperature of the cooling water increases. The net heat rejected is given by Q_re = h₄ – h_1.

At constant pressure, the enthalpy change equals the energy transferred from the environment through heating:

Isobaric process (Vdp = 0):

dh = dQ → Q = h₂ – h₁

At constant entropy, i.e. in isentropic process, the enthalpy change equals the flow process work done on or by the system:

Isentropic process (dQ = 0):

dh = Vdp → W = h₂ – h₁

W_net = W_T - W_P = (h₄ – h₃)-(h₂ – h₁)

Q_add = h₃ – h₂

`eta` = W_net/Q_add

`eta`= (Q_add - Q_re)/Q_re

`eta` = 1 - (Q_re/Q_add)

Required Inputs:

Turbine Inlet temperature = 400dC

Turbine Inlet Pressure = 30 bars

Turbine Outlet Pressure = 0.05 bars = Condenser Inlet Pressure

Matlab Script:

close all
clear all
clc

% Print Process Description
fprintf('1-2 Isentropic compression\n')
fprintf('2-3 Isobaric heat addition\n')
fprintf('3-4 Isentropic expansion\n')
fprintf('4-1 Isobaric heat rejection\n')

%step 1: Enter Inputs
p1 = input('Enter the pressure at the turbine inlet(in bar):')
T1 = input('Enter the Temperature at the turbine inlet(in deg):')
p2 = input('Enter the pressure at the Condenser (in bar):')

%step 2: Calculate each state point values.
% state point 1
h1 = XSteam('h_pt',p1,T1)
s1 = XSteam('s_pt',p1,T1)

%step 3: state point 2
s2 = s1
sf2 = XSteam('sL_p',p2)
sg2 = XSteam('sV_p',p2)
hf2 = XSteam('hL_p',p2)
hg2 = XSteam('hV_p',p2)
x2 = (s2-sf2)/(sg2-sf2)
s2 = sg2
h2 = hg2 %hf2+(x2*(hg2-hf2))
T2 = XSteam('T_hs',h2,s2)

%step 4: state point 3
T3 = T2
p3 = p2
s3 = XSteam('sL_T',T3)
h3 = XSteam('hL_T',T3)
v3 = XSteam('vL_p',p2)

%step 5: state point 4
p4 = p1
s4 = s3
Cp4 = XSteam('Cp_ps',p4,s4)
Cv4 = XSteam('Cv_ps',p4,s4)
gamma = Cp4/Cv4
r = (gamma - 1)/gamma
T4 = T3*(p4/p3)^r

%step 6: Calculate other output values.
%Turbine Work
Wt = h1-h2
Wp = v3*(p1-p2)*100  %Pump Work
h4 = Wp + h3
W_net = Wt-Wp  % Net Work unit in kj/kg
Qs = h1 - h4 %heat supplied
N_per = (W_net/Qs)*100 %Efficiency
steam_rate = 3600/(Wt - Wp)  %Steam rate per second

%step 7: state pont 5 and 6 which help in plotting
T5 = XSteam('Tsat_p',p1)
T6 = T5
s5 = XSteam('sL_p',p1)
s6 = XSteam('sV_p',p1)
h5 = XSteam('hL_p',p1)
h6 = XSteam('hV_p',p1)

%step 8: print all state point values to shaow.
fprintf('\n--------------RESULT-------------\n')
fprintf('At State point 1:\nP1 is : %4.2f Bar\nT1 is : %4.2f Degree Celcius\nH1 is : %4.2f kJ/kg\nS1 is : %4.2f kJ/kgk\n\n',p1,T1,h1,s1)
fprintf('At State point 2:\nP2 is : %4.2f Bar\nT2 is : %4.2f Degree Celcius\nH2 is : %4.2f kJ/kg\nS2 is : %4.2f kJ/kgk\n\n',p1,T2,h2,s2)
fprintf('At State point 3:\nP3 is : %4.2f Bar\nT3 is : %4.2f Degree Celcius\nH3 is : %4.2f kJ/kg\nS3 is : %4.2f kJ/kgk\n\n',p1,T3,h3,s3)
fprintf('At State point 4:\nP4 is : %4.2f Bar\nT4 is : %4.2f Degree Celcius\nH4 is : %4.2f kJ/kg\nS4 is : %4.2f kJ/kgk\n\n',p1,T4,h4,s4)

%step 9: plot T-S Daigram
subplot(1,2,1)
T_rankine = linspace(T5,T1,5000);
for i = 1:length(T_rankine);
    s_rankine(i) = XSteam('s_pT',p1,T_rankine(i));
end
plot([s1 s2 s3 s4 s5 s6 s_rankine],[T1 T2 T3 T4 T5 T6 T_rankine],'linewidth',1.5)
hold on
T_curve = linspace(1,400,500);
for i=1:length(T_curve)
    s_L(i) = XSteam('sL_T',T_curve(i));
    s_V(i) = XSteam('sV_T',T_curve(i));
    h_L(i) = XSteam('hL_T',T_curve(i));
    h_V(i) = XSteam('hV_T',T_curve(i));
end

plot(s_L,T_curve)
plot(s_V,T_curve)
text(s1,T1,'1','fontsize',15)
text(s2,T2,'2','fontsize',15)
text(s3,T3,'3','fontsize',15,'horizontalalignment','right','verticalalignment','top')
text(s4,T4,'4','fontsize',15,'horizontalalignment','right','verticalalignment','bottom')
text(s5,T5,'5','fontsize',15)
text(s6,T6,'6','fontsize',15)
xlabel('Entropy in kJ/KgK')
ylabel('Temperature in Degree Celcius')

%step 10: plot H-S diagram
subplot(1,2,2)
plot(s_V,h_V,'color','r')
plot(s_L,h_L,'color','r')
text(s1,T1,'1','fontsize',15)
text(s2,T2,'2','fontsize',15)
text(s3,T3,'3','fontsize',15,'horizontalalignment','right','verticalalignment','top')
text(s4,T4,'4','fontsize',15,'horizontalalignment','right','verticalalignment','bottom')
xlabel('Entropy in kJ/KgK')
ylabel('Enthalpy in kJ/kg')

Results: