# -*- 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()