# -*- coding: utf-8 -*-
"""
Created on Wed Feb 26 11:11:00 2025

@author: Ronelly De Souza
"""

#Ideal Rankine cycle with Reheat
#Objectives:
#calculate thermodynamic properties in each point of the cycle
#calculate the power generated by the turbine
#calculate the power consumed by the pump
#calculate the heat supplied at the vapor generator
#calculate the heat rejected at the condenser
#calculate the thermal efficiency of the cycle
#plot the water T-s diagram with the modelled Rankine cycle

import numpy as np
import matplotlib.pyplot as plt
from CoolProp.CoolProp import *

#input data
P_net = 1E5 #kW

c_Pa = 1E5 #factor to convert from bar to Pa

T_sat = np.linspace(273.15,647,100) #limited by the temperature range defined within CoolProp database
s_liq = [PropsSI("S","T",T,"Q",0,"Water") for T in T_sat] #calculates the entropy at the saturated liquid point for each temperature in T_sat
s_vap = [PropsSI("S","T",T,"Q",1,"Water") for T in T_sat] #calculates the entropy at the saturated vapor point for each temperature in T_sat

print("Ideal Rankine Cycle model with Reheat")

#Point 1: turbine inlet
P1 = 80 #bar
T1 = 480+273 #K
h1 = PropsSI("H", "P", P1*c_Pa, "T", T1, "Water")
s1 = PropsSI("S", "P", P1*c_Pa, "T", T1, "Water")

print()
print("Properties at point 1")
print(f"P1 = {'%.2f'%(P1)} bar")
print(f"T1 = {'%.2f'%(T1-273.15)} °C")
print(f"h1 = {'%.2f'%(h1/1000)} kJ/kg")
print(f"s1 = {'%.2f'%(s1/1000)} kJ/kg.K")

#Point 2: first stage turbine outlet
P2 = 7 #bar
s2 = s1 #entropy for isentropic expansion
x2 = PropsSI("Q", "P", P2*c_Pa, "S", s2, "Water")
h2 = PropsSI("H", "P", P2*c_Pa, "S", s2, "Water")
T2 = PropsSI("T", "P", P2*c_Pa, "S", s2, "Water")

print()
print("Properties at point 2")
print(f"P2 = {'%.2f'%(P2)} bar")
print(f"T2 = {'%.2f'%(T2-273.15)} °C")
print(f"h2 = {'%.2f'%(h2/1000)} kJ/kg")
print(f"s2 = {'%.2f'%(s2/1000)} kJ/kg.K")
print(f"x2 = {'%.2f'%(x2)}")

#Point 3: second stage turbine inlet
P3 = P2 #bar
T3 = 440+273 #K
h3 = PropsSI("H", "P", P3*c_Pa, "T", T3, "Water")
s3 = PropsSI("S", "P", P3*c_Pa, "T", T3, "Water")

print()
print("Properties at point 3")
print(f"P3 = {'%.2f'%(P3)} bar")
print(f"T3 = {'%.2f'%(T3-273.15)} °C")
print(f"h3 = {'%.2f'%(h3/1000)} kJ/kg")
print(f"s3 = {'%.2f'%(s3/1000)} kJ/kg.K")

#Point 4: second stage turbine outlet or condenser inlet
P4 = 0.08 #bar
s4 = s3 #entropy for isentropic expansion
h4 = PropsSI("H", "P", P4*c_Pa, "S", s4, "Water")
x4 = PropsSI("Q", "P", P4*c_Pa, "H", h4, "Water")
T4 = PropsSI("T", "P", P4*c_Pa, "S", s4, "Water")

print()
print("Properties at point 4")
print(f"P4 = {'%.2f'%(P4)} bar")
print(f"s4 = {'%.2f'%(s4/1000)} kJ/kg.K")
print(f"h4 = {'%.2f'%(h4/1000)} kJ/kg")
print(f"x4 = {'%.2f'%(x4)}")
print(f"T4 = {'%.2f'%(T4-273.15)} °C")

#point 5: condenser outlet
P5 = P4
x5 = 0 #quality
h5 = PropsSI("H", "P", P5*c_Pa, "Q", x5, "Water")
s5 = PropsSI("S", "P", P5*c_Pa, "Q", x5, "Water")
T5 = PropsSI("T", "P", P5*c_Pa, "Q", x5, "Water")

print()
print("Properties at point 5")
print(f"P5 = {'%.2f'%(P5)} bar")
print(f"x5 = {'%.2f'%(x5)}")
print(f"h5 = {'%.2f'%(h5/1000)} kJ/kg")
print(f"s5 = {'%.2f'%(s5/1000)} kJ/kg.K")
print(f"T5 = {'%.2f'%(T5-273.15)} °C")

#point 6: pump outlet
P6 = P1
s6 = s5
h6 = PropsSI("H", "P", P6*c_Pa, "S", s6, "Water")
T6 = PropsSI("T", "P", P6*c_Pa, "S", s6, "Water")

print()
print("Properties at point 6")
print(f"P5 = {'%.2f'%(P6)} bar")
print(f"s6 = {'%.2f'%(s6/1000)} kJ/kg.K")
print(f"h6 = {'%.2f'%(h6/1000)} kJ/kg")
print(f"T6 = {'%.2f'%(T6-273.15)} °C")

#point 7: beginning of evaporation
P7 = P1
x7 = 0
h7 = PropsSI("H", "P", P7*c_Pa, "Q", x7, "Water")
s7 = PropsSI("S", "P", P7*c_Pa, "Q", x7, "Water")
T7 = PropsSI("T", "P", P7*c_Pa, "Q", x7, "Water")

print()
print("Properties at point 7")
print(f"P7 = {'%.2f'%(P7)} bar")
print(f"s7 = {'%.2f'%(s7/1000)} kJ/kg.K")
print(f"h7 = {'%.2f'%(h7/1000)} kJ/kg")
print(f"T7 = {'%.2f'%(T7-273.15)} °C")

#point 8: ending of evaporation
P8 = P1
x8 = 1
h8 = PropsSI("H", "P", P8*c_Pa, "Q", x8, "Water")
s8 = PropsSI("S", "P", P8*c_Pa, "Q", x8, "Water")
T8 = PropsSI("T", "P", P8*c_Pa, "Q", x8, "Water")

print()
print("Properties at point 8")
print(f"P8 = {'%.2f'%(P8)} bar")
print(f"s8 = {'%.2f'%(s8/1000)} kJ/kg.K")
print(f"h8 = {'%.2f'%(h8/1000)} kJ/kg")
print(f"T8 = {'%.2f'%(T8-273.15)} °C")


#calculate the power generated by the turbine
m_vap = P_net/((h1-h2)+(h3-h4)-(h6-h5)) #kg/s
P_t1 = m_vap*(h1-h2) #kW
P_t2 = m_vap*(h3-h4) #kW

#calculate the power consumed by the pump
P_p = m_vap*(h6-h5)

#calculate the heat supplied at the vapor generator
heat_in = m_vap*((h1-h6)+(h3-h2))

#calculate the heat rejected at the condenser
heat_out = m_vap*(h4-h5)

#calculate the thermal efficiency of the cycle
eta_th = P_net/heat_in

print()
print("Power generated by the turbine stage 1:")
print(f"P_t1 = {'%.2f'%(P_t1)} kW")
print()
print("Power generated by the turbine stage 2:")
print(f"P_t2 = {'%.2f'%(P_t2)} kW")
print()
print("Power consumed by the pump:")
print(f"P_p = {'%.2f'%(P_p)} kW")
print()
print("Net power from the cycle:")
print(f"P_net = {'%.2f'%(P_net)} kW")
print()
print("Heat supplied at the vapor generator:")
print(f"heat_in = {'%.2f'%(heat_in)} kW")
print()
print("Heat rejected at the condenser:")
print(f"heat_out = {'%.2f'%(heat_out)} kW")
print()
print("Thermal efficiency of the cycle")
print(f"eta_th = {'%.2f'%(eta_th*100)}%")

plt.plot(s_liq, T_sat, "b-", label="Saturation Liquid Line")  #liquid saturation line
plt.plot(s_vap, T_sat, "r-", label="Saturation Vapor Line")   #vapor saturation line

plt.xlabel("Entropy (J/kg·K)")
plt.ylabel("Temperature (K)")
plt.title("Water T-s Diagram with Rankine Cycle")

s_cycle = [s1, s2, s3, s4, s5, s6, s7, s8, s1] #entropy for each point
T_cycle = [T1, T2, T3, T4, T5, T6, T7, T8, T1] #temperature for each point
plt.plot(s_cycle, T_cycle, "ko--", label="Rankine Cycle with Reheat")  #plot the cycle points

plt.legend()

#"for" and "if" loops to show the number of each thermodynamic state in an appropriate position
for i in range(len(s_cycle)): #iterates each value of the s_cycle list. Thermodynamic state "1" starts at python index "0"
    if i == 0:
        x_offset = 100 #value to move the text horizontally with respect to the s_cycle[i] value
        y_offset = -10 #value to move the text vertically with respect to the s_cycle[i] value
        plt.text(s_cycle[i]+x_offset, T_cycle[i]+y_offset, i + 1)
    elif i == 4:
        x_offset = 100
        y_offset = -20
        plt.text(s_cycle[i]+x_offset, T_cycle[i]+y_offset, f"point {i+1}")
    elif i == 8:
        plt.text(s_cycle[i]+x_offset, T_cycle[i]+y_offset, "") #this would be the point 9, which does not need to be shown
    else:
        x_offset = -300
        y_offset = 5
        plt.text(s_cycle[i]+x_offset, T_cycle[i]+y_offset, i + 1)

plt.show()