program rankine_regenerative implicit none integer :: fwh_type, n_extractions, i, iostat_val double precision :: P_boiler, T_boiler, P_cond, P_ext1, P_ext2 double precision :: eta_turbine, eta_pump double precision :: h1, s1, h2a, h2s, h3a, h3s, h4, h5, h6a, h6s double precision :: h_f_cond, h_fg_cond, s_f_cond, s_fg_cond, T_sat_cond double precision :: h_f_ext1, h_fg_ext1, s_f_ext1, s_fg_ext1, T_sat_ext1 double precision :: h_f_ext2, h_fg_ext2, s_f_ext2, s_fg_ext2, T_sat_ext2 double precision :: h_f_boil, s_f_boil, T_sat_boil, v_f_cond, v_f_ext1 double precision :: x2, x3, y1, y2, Wt, Wp, Wnet, Qin, Qout double precision :: eta_th, bwr, m_dot, P_output, SSC double precision :: h_fw, h_mix, v_f_ext2 double precision :: P_i, eta_i, y_i double precision, parameter :: Ru = 8.314462d0 character(len=40) :: fwh_name read(*,*,iostat=iostat_val) P_boiler if (iostat_val /= 0) then write(*,*) 'ERROR: Invalid boiler pressure input.' stop end if read(*,*,iostat=iostat_val) T_boiler read(*,*,iostat=iostat_val) P_cond read(*,*,iostat=iostat_val) n_extractions read(*,*,iostat=iostat_val) P_ext1 read(*,*,iostat=iostat_val) P_ext2 read(*,*,iostat=iostat_val) fwh_type read(*,*,iostat=iostat_val) eta_turbine read(*,*,iostat=iostat_val) eta_pump read(*,*,iostat=iostat_val) m_dot if (iostat_val /= 0) then write(*,*) 'ERROR: Failed to read all Rankine inputs.' stop end if if (P_boiler <= 0.0d0 .or. P_cond <= 0.0d0) then write(*,*) 'ERROR: Pressures must be positive.' stop end if if (P_boiler <= P_cond) then write(*,*) 'ERROR: Boiler pressure must exceed condenser pressure.' stop end if if (eta_turbine <= 0.0d0 .or. eta_turbine > 1.0d0) eta_turbine = 0.88d0 if (eta_pump <= 0.0d0 .or. eta_pump > 1.0d0) eta_pump = 0.85d0 if (n_extractions < 1) n_extractions = 1 if (n_extractions > 2) n_extractions = 2 if (P_ext1 <= P_cond .or. P_ext1 >= P_boiler) P_ext1 = sqrt(P_boiler*P_cond) if (n_extractions == 2) then if (P_ext2 <= P_cond .or. P_ext2 >= P_ext1) & P_ext2 = sqrt(P_cond*P_ext1) end if if (m_dot <= 0.0d0) m_dot = 1.0d0 if (fwh_type == 1) then fwh_name = 'Open feedwater heater (direct contact)' else fwh_name = 'Closed feedwater heater (shell-tube)' fwh_type = 2 end if ! Compute saturation properties at key pressures call sat_props(P_cond, T_sat_cond, h_f_cond, h_fg_cond, s_f_cond, s_fg_cond, v_f_cond) call sat_props(P_ext1, T_sat_ext1, h_f_ext1, h_fg_ext1, s_f_ext1, s_fg_ext1, v_f_ext1) if (n_extractions == 2) then call sat_props(P_ext2, T_sat_ext2, h_f_ext2, h_fg_ext2, s_f_ext2, s_fg_ext2, v_f_ext2) end if call sat_props(P_boiler, T_sat_boil, h_f_boil, h_fg_cond, s_f_boil, s_fg_cond, v_f_cond) ! State 1: Turbine inlet (superheated steam at P_boiler, T_boiler) call superheat_props(P_boiler, T_boiler, T_sat_boil, h_f_boil, s_f_boil, h1, s1) ! State 2: After first extraction (isentropic expansion to P_ext1) call sat_props(P_ext1, T_sat_ext1, h_f_ext1, h_fg_ext1, s_f_ext1, s_fg_ext1, v_f_ext1) x2 = (s1 - s_f_ext1) / max(s_fg_ext1, 1.0d-10) if (x2 > 1.0d0) then ! Still superheated at extraction call superheat_from_entropy(P_ext1, s1, T_sat_ext1, h_f_ext1, s_f_ext1, s_fg_ext1, h2s) else h2s = h_f_ext1 + x2 * h_fg_ext1 end if h2a = h1 - eta_turbine * (h1 - h2s) ! State 3: Turbine exit (isentropic expansion to P_cond) call sat_props(P_cond, T_sat_cond, h_f_cond, h_fg_cond, s_f_cond, s_fg_cond, v_f_cond) x3 = (s1 - s_f_cond) / max(s_fg_cond, 1.0d-10) if (x3 > 1.0d0) then call superheat_from_entropy(P_cond, s1, T_sat_cond, h_f_cond, s_f_cond, s_fg_cond, h3s) else h3s = h_f_cond + x3 * h_fg_cond end if h3a = h1 - eta_turbine * (h1 - h3s) ! State 4: Condenser exit (saturated liquid at P_cond) h4 = h_f_cond ! Pump work and extraction fraction calculation if (n_extractions == 1) then if (fwh_type == 1) then ! OPEN FWH: 1 extraction ! Pump 1: P_cond → P_ext1 call sat_props(P_cond, T_sat_cond, h_f_cond, h_fg_cond, s_f_cond, s_fg_cond, v_f_cond) h5 = h4 + v_f_cond*(P_ext1 - P_cond)/eta_pump ! after pump 1 ! FWH energy balance: y*h2a + (1-y)*h5 = h_f_ext1 h_fw = h_f_ext1 y1 = (h_fw - h5) / max(h2a - h5, 1.0d-10) if (y1 < 0.0d0) y1 = 0.0d0 if (y1 > 1.0d0) y1 = 1.0d0 ! Pump 2: P_ext1 → P_boiler h6a = h_fw + v_f_ext1*(P_boiler - P_ext1)/eta_pump ! Turbine work per unit mass at boiler Wt = (h1 - h2a) + (1.0d0 - y1)*(h2a - h3a) ! Pump work Wp = (1.0d0 - y1)*(h5 - h4) + (h6a - h_fw) Qin = h1 - h6a else ! CLOSED FWH: 1 extraction h5 = h4 + v_f_cond*(P_boiler - P_cond)/eta_pump ! Closed FWH: extracted steam heats feedwater ! y*h2a + h5 = y*h_f_ext1 + h_fw (drain to trap) ! Assume drain cascaded back: h_fw = h5 + y*(h2a - h_f_ext1) y1 = (h_f_ext1 - h5) / max((h2a - h_f_ext1) + (h_f_ext1 - h5), 1.0d-10) if (y1 < 0.0d0) y1 = 0.0d0 if (y1 > 1.0d0) y1 = 1.0d0 h6a = h5 + y1*(h2a - h_f_ext1) Wt = (h1 - h2a) + (1.0d0 - y1)*(h2a - h3a) Wp = (h5 - h4) Qin = h1 - h6a end if y2 = 0.0d0 else ! 2 extractions — open FWH simplified call sat_props(P_ext2, T_sat_ext2, h_f_ext2, h_fg_ext2, s_f_ext2, s_fg_ext2, v_f_ext2) ! Intermediate extraction state ! For simplicity: second extraction between ext1 and condenser ! Recompute h at P_ext2 from isentropic expansion x2 = (s1 - s_f_ext2) / max(s_fg_ext2, 1.0d-10) if (x2 > 1.0d0) then call superheat_from_entropy(P_ext2, s1, T_sat_ext2, h_f_ext2, & s_f_ext2, s_fg_ext2, h_mix) else h_mix = h_f_ext2 + x2 * h_fg_ext2 end if h_mix = h1 - eta_turbine*(h1 - h_mix) ! actual h at ext2 ! Pump 1: P_cond → P_ext2 h5 = h4 + v_f_cond*(P_ext2 - P_cond)/eta_pump ! Open FWH 1 at P_ext2 y2 = (h_f_ext2 - h5) / max(h_mix - h5, 1.0d-10) if (y2 < 0.0d0) y2 = 0.0d0 if (y2 > 0.5d0) y2 = 0.5d0 ! Pump 2: P_ext2 → P_ext1 h_fw = h_f_ext2 + v_f_ext2*(P_ext1 - P_ext2)/eta_pump ! Open FWH 2 at P_ext1 y1 = (h_f_ext1 - h_fw) / max(h2a - h_fw, 1.0d-10) if (y1 < 0.0d0) y1 = 0.0d0 if (y1 > 0.5d0) y1 = 0.5d0 ! Pump 3: P_ext1 → P_boiler h6a = h_f_ext1 + v_f_ext1*(P_boiler - P_ext1)/eta_pump Wt = (h1 - h2a) + (1.0d0-y1)*(h2a - h_mix) + (1.0d0-y1-y2)*(h_mix - h3a) Wp = (1.0d0-y1-y2)*(h5-h4) + (1.0d0-y1)*(h_fw-h_f_ext2) + (h6a-h_f_ext1) Qin = h1 - h6a end if Wnet = Wt - Wp Qout = (1.0d0 - y1 - y2) * (h3a - h4) eta_th = Wnet / max(Qin, 1.0d-10) bwr = Wp / max(Wt, 1.0d-10) P_output = m_dot * Wnet / 1000.0d0 ! kW SSC = 3600.0d0 / max(Wnet, 1.0d-10) ! kg/kWh write(*,'(A)') '============================================================' write(*,'(A)') ' RANKINE CYCLE — REGENERATIVE ENGINE' write(*,'(A)') '============================================================' write(*,*) write(*,'(A)') '--- INPUTS --------------------------------------------------' write(*,'(A,ES12.4,A)') ' Boiler Pressure = ', P_boiler, ' Pa' write(*,'(A,F12.2,A)') ' Boiler Temperature = ', T_boiler, ' K' write(*,'(A,ES12.4,A)') ' Condenser Pressure = ', P_cond, ' Pa' write(*,'(A,I8)') ' Number of Extractions = ', n_extractions write(*,'(A,ES12.4,A)') ' Extraction 1 Pressure = ', P_ext1, ' Pa' if (n_extractions == 2) & write(*,'(A,ES12.4,A)') ' Extraction 2 Pressure = ', P_ext2, ' Pa' write(*,'(A,A)') ' FWH Type = ', trim(fwh_name) write(*,'(A,F10.4)') ' Turbine Isentropic Eff = ', eta_turbine write(*,'(A,F10.4)') ' Pump Isentropic Eff = ', eta_pump write(*,'(A,ES12.4,A)') ' Mass Flow Rate = ', m_dot, ' kg/s' write(*,*) write(*,'(A)') '--- SATURATION DATA -----------------------------------------' write(*,'(A,F12.2,A)') ' T_sat (boiler) = ', T_sat_boil, ' K' write(*,'(A,F12.2,A)') ' T_sat (ext1) = ', T_sat_ext1, ' K' if (n_extractions == 2) & write(*,'(A,F12.2,A)') ' T_sat (ext2) = ', T_sat_ext2, ' K' write(*,'(A,F12.2,A)') ' T_sat (condenser) = ', T_sat_cond, ' K' write(*,*) write(*,'(A)') '--- STATE POINTS (per kg at boiler) -------------------------' write(*,'(A)') ' Point h[kJ/kg] s[kJ/kgK] T[K] P[Pa]' write(*,'(A)') ' ----------------------------------------------------------------' write(*,'(A,F12.2,2X,F10.4,2X,F12.2,2X,ES12.4)') ' 1 turb in ', h1/1000, s1/1000, T_boiler, P_boiler write(*,'(A,F12.2,2X,A,2X,F12.2,2X,ES12.4)') ' 2 ext1 ', h2a/1000, ' --- ', T_sat_ext1, P_ext1 write(*,'(A,F12.2,2X,A,2X,F12.2,2X,ES12.4)') ' 3 turb out ', h3a/1000, ' --- ', T_sat_cond, P_cond write(*,'(A,F12.2,2X,A,2X,F12.2,2X,ES12.4)') ' 4 cond out ', h4/1000, ' --- ', T_sat_cond, P_cond write(*,'(A,F12.2,2X,A,2X,A,2X,ES12.4)') ' 6 boil in ', h6a/1000, ' --- ', ' --- ', P_boiler write(*,*) write(*,'(A)') '--- EXTRACTION FRACTIONS ------------------------------------' write(*,'(A,F10.5)') ' y1 (extraction 1) = ', y1 if (n_extractions == 2) & write(*,'(A,F10.5)') ' y2 (extraction 2) = ', y2 write(*,'(A,F10.5)') ' Condenser fraction = ', 1.0d0-y1-y2 write(*,*) write(*,'(A)') '--- ENERGY PER UNIT MASS ------------------------------------' write(*,'(A,F12.2,A)') ' Turbine Work Wt = ', Wt/1000, ' kJ/kg' write(*,'(A,F12.2,A)') ' Pump Work Wp = ', Wp/1000, ' kJ/kg' write(*,'(A,F12.2,A)') ' Net Work Wnet = ', Wnet/1000, ' kJ/kg' write(*,'(A,F12.2,A)') ' Heat Input Qin = ', Qin/1000, ' kJ/kg' write(*,'(A,F12.2,A)') ' Heat Rejected Qout = ', Qout/1000, ' kJ/kg' write(*,*) write(*,'(A)') '--- CYCLE PERFORMANCE ---------------------------------------' write(*,'(A,F10.4)') ' Thermal Efficiency = ', eta_th write(*,'(A,F10.2,A)') ' Thermal Efficiency = ', eta_th*100.0d0, ' percent' write(*,'(A,F10.4)') ' Back-Work Ratio = ', bwr write(*,'(A,F12.2,A)') ' Power Output = ', P_output, ' kW' write(*,'(A,F10.4,A)') ' Specific Steam Consumption= ', SSC, ' kg/kWh' write(*,*) ! Efficiency vs extraction pressure sweep write(*,'(A)') '--- EFFICIENCY VS EXTRACTION PRESSURE SWEEP -----------------' write(*,'(A)') ' P_ext1[kPa] y1 eta_th Wnet[kJ/kg]' write(*,'(A)') ' -----------------------------------------------------------' do i = 1, 40 P_i = P_cond + (P_boiler - P_cond) * dble(i) / 41.0d0 call compute_single_extraction(P_boiler, T_boiler, P_cond, P_i, & eta_turbine, eta_pump, h1, s1, h_f_cond, v_f_cond, eta_i, y_i) write(*,'(F12.2,2X,F10.5,2X,F10.5,2X,F12.2)') P_i/1000.0d0, y_i, eta_i, & eta_i * (h1 - h_f_cond - v_f_cond*(P_boiler-P_cond)/eta_pump) / 1000.0d0 end do write(*,*) write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------' write(*,'(A)') ' Simplified water saturation correlations (curve fits).' write(*,'(A)') ' Isentropic turbine: h_out = h_in - eta_t*(h_in - h_out_s).' write(*,'(A)') ' Pump work: w_p = v_f*(P_out - P_in)/eta_p.' write(*,'(A)') ' Open FWH: y*h_ext + (1-y)*h_fw = h_f_sat.' contains subroutine sat_props(P_pa, Tsat, hf, hfg, sf, sfg, vf) implicit none double precision, intent(in) :: P_pa double precision, intent(out) :: Tsat, hf, hfg, sf, sfg, vf double precision :: P_MPa, logP P_MPa = P_pa / 1.0d6 if (P_MPa < 0.001d0) P_MPa = 0.001d0 logP = log(P_MPa) ! Curve-fit approximations for water (valid ~5 kPa to 22 MPa) Tsat = 373.15d0 + 42.0d0*logP - 0.8d0*logP**2 ! K if (Tsat < 290.0d0) Tsat = 290.0d0 + 15.0d0*P_MPa if (Tsat > 647.0d0) Tsat = 647.0d0 hf = (417.0d0 + 390.0d0*logP + 35.0d0*logP**2) * 1000.0d0 ! J/kg if (hf < 100.0d3) hf = 100.0d3 + 200.0d3*P_MPa if (hf > 2100.0d3) hf = 2100.0d3 hfg = (2258.0d0 - 180.0d0*logP - 40.0d0*logP**2) * 1000.0d0 if (hfg < 200.0d3) hfg = 200.0d3 if (hfg > 2500.0d3) hfg = 2500.0d3 sf = (1.303d0 + 0.82d0*logP + 0.06d0*logP**2) * 1000.0d0 ! J/(kg K) if (sf < 0.3d3) sf = 0.3d3 + 0.5d3*P_MPa sfg = hfg / Tsat vf = 0.001d0 * (1.0d0 + 0.0002d0*P_MPa) ! m3/kg end subroutine sat_props subroutine superheat_props(P_pa, T, Tsat, hf_boil, sf_boil, h_out, s_out) implicit none double precision, intent(in) :: P_pa, T, Tsat, hf_boil, sf_boil double precision, intent(out) :: h_out, s_out double precision :: P_MPa, hfg_b, sfg_b, hg, sg, cp_steam, dT P_MPa = P_pa / 1.0d6 hfg_b = (2258.0d0 - 180.0d0*log(P_MPa) - 40.0d0*log(P_MPa)**2)*1000.0d0 if (hfg_b < 200.0d3) hfg_b = 200.0d3 hg = hf_boil + hfg_b sfg_b = hfg_b / Tsat sg = sf_boil + sfg_b cp_steam = 2100.0d0 ! J/(kg K) approximate dT = T - Tsat if (dT < 0.0d0) dT = 0.0d0 h_out = hg + cp_steam * dT s_out = sg + cp_steam * log(max(T/Tsat, 1.0d0)) end subroutine superheat_props subroutine superheat_from_entropy(P_pa, s_target, Tsat, hf, sf, sfg, h_out) implicit none double precision, intent(in) :: P_pa, s_target, Tsat, hf, sf, sfg double precision, intent(out) :: h_out double precision :: hfg, hg, sg, cp_steam, T_est, dT hfg = sfg * Tsat hg = hf + hfg sg = sf + sfg cp_steam = 2100.0d0 ! s = sg + cp ln(T/Tsat) => T = Tsat * exp((s-sg)/cp) T_est = Tsat * exp((s_target - sg)/cp_steam) dT = T_est - Tsat if (dT < 0.0d0) dT = 0.0d0 h_out = hg + cp_steam * dT end subroutine superheat_from_entropy subroutine compute_single_extraction(Pb, Tb, Pc, Pe, etat, etap, & h1in, s1in, hfc, vfc, eta_out, y_out) implicit none double precision, intent(in) :: Pb, Tb, Pc, Pe, etat, etap, h1in, s1in, hfc, vfc double precision, intent(out) :: eta_out, y_out double precision :: Tse, hfe, hfge, sfe, sfge, vfe, xe, h2si, h2ai double precision :: h3si, h3ai, x3i, h5i, h6i, Wti, Wpi, Qini double precision :: Tsc, hfc_local, hfgc, sfc, sfgc, vfc2 call sat_props(Pe, Tse, hfe, hfge, sfe, sfge, vfe) call sat_props(Pc, Tsc, hfc_local, hfgc, sfc, sfgc, vfc2) xe = (s1in - sfe)/max(sfge,1.0d-10) if (xe > 1.0d0) then call superheat_from_entropy(Pe, s1in, Tse, hfe, sfe, sfge, h2si) else h2si = hfe + xe*hfge end if h2ai = h1in - etat*(h1in - h2si) x3i = (s1in - sfc)/max(sfgc,1.0d-10) if (x3i > 1.0d0) then call superheat_from_entropy(Pc, s1in, Tsc, hfc_local, sfc, sfgc, h3si) else h3si = hfc_local + x3i*hfgc end if h3ai = h1in - etat*(h1in - h3si) h5i = hfc_local + vfc2*(Pe-Pc)/etap y_out = (hfe - h5i)/max(h2ai - h5i, 1.0d-10) if (y_out < 0.0d0) y_out = 0.0d0 if (y_out > 1.0d0) y_out = 1.0d0 h6i = hfe + vfe*(Pb-Pe)/etap Wti = (h1in-h2ai)+(1.0d0-y_out)*(h2ai-h3ai) Wpi = (1.0d0-y_out)*(h5i-hfc_local)+(h6i-hfe) Qini = h1in - h6i eta_out = (Wti-Wpi)/max(Qini,1.0d-10) end subroutine compute_single_extraction end program rankine_regenerative