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

@author: Ronelly De Souza
"""

#Rankine cycle with irreversibilities
#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 *


# Generate saturation lines (liquid and vapor)
T_sat = np.linspace(273.06, 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("Rankine Cycle model with irreversibilities")

#input data
P_net = 1E5 #kW
eta_t = 0.85 #turbine isentropic efficiency
eta_p = 0.85 #turbine isentropic efficiency

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

#Point 1: turbine inlet
P1 = 80 #bar
x1 = 1 #quality
T1 = PropsSI("T", "P", P1*c_Pa, "Q", x1, "Water")
h1 = PropsSI("H", "P", P1*c_Pa, "Q", x1, "Water")
s1 = PropsSI("S", "P", P1*c_Pa, "Q", x1, "Water")

print()
print("Properties at point 1")
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")
print(f"x1 = {x1}")

#Point 2: turbine outlet
P2 = 0.08 #bar
s2_is = s1 #entropy for an isentropic expansion
h2_is = PropsSI("H", "P", P2*c_Pa, "S", s2_is, "Water")
x2_is = PropsSI("Q", "P", P2*c_Pa, "S", s2_is, "Water")
h2 = h1 - ((h1 - h2_is)*eta_t)
x2 = PropsSI("Q", "P", P2*c_Pa, "H", h2, "Water")
s2 = PropsSI("S", "P", P2*c_Pa, "H", h2, "Water")
T2 = PropsSI("T", "P", P2*c_Pa, "S", s2, "Water")

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

#Point 3: pump inlet
P3 = P2 #bar
x3 = 0 #quality
T3 = PropsSI("T", "P", P3*c_Pa, "Q", x3, "Water")
h3 = PropsSI("H", "P", P3*c_Pa, "Q", x3, "Water")
s3 = PropsSI("S", "P", P3*c_Pa, "Q", x3, "Water")

print()
print("Properties at point 3")
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")
print(f"x3 = {x3}")

#Point 4: pump outlet
P4 = P1 #bar
s4_is = s3 #entropy for an isentropic expansion
h4_is = PropsSI("H", "P", P4*c_Pa, "S", s4_is, "Water")
h4 = h3 + ((h4_is - h3)/eta_t)
s4 = PropsSI("S", "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"T4 = {'%.2f'%(T4-273.15)} °C")
print(f"h4 isentropic = {'%.2f'%(h4_is/1000)} kJ/kg")
print(f"h4 = {'%.2f'%(h4/1000)} kJ/kg")
print(f"s4 isentropic = {'%.2f'%(s4_is/1000)} kJ/kg")
print(f"s4 = {'%.2f'%(s4/1000)} kJ/kg.K")

#Point 5: evaporator inlet
P5 = P1
x5 = 0 #quality
T5 = T1
h5 = PropsSI("H", "P", P5*c_Pa, "Q", x5, "Water")
s5 = PropsSI("S", "P", P5*c_Pa, "Q", x5, "Water")

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

#calculate the power generated by the turbine
m_vap = P_net/((h1-h2)-(h4-h3)) #kg/s
P_t = m_vap*(h1-h2) #kW

#calculate the power consumed by the pump
P_p = m_vap*(h4-h3)

#calculate the heat supplied at the vapor generator
heat_in = m_vap*(h1-h4)

#calculate the heat rejected at the condenser
heat_out = m_vap*(h3-h2)

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

print()
print("Power generated by the turbine:")
print(f"P_t = {'%.2f'%(P_t)} 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)}%")

#Plot T-s Diagram
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")

#plot the thermodynamic states of the modelled Rankine Cycle
s_cycle = [s1, s2, s3, s4, s5, s1] #entropy for each point
T_cycle = [T1, T2, T3, T4, T5, T1] #temperature for each point
plt.plot(s_cycle, T_cycle, "ko--", label="Rankine Cycle")  #plot the cycle points

plt.legend(loc="upper right", fontsize=9)

#"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 = 200 #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 == 2:
        x_offset = 100
        y_offset = -20
        plt.text(s_cycle[i]+x_offset, T_cycle[i]+y_offset, i + 1)
    elif i == 5:
        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()