program combined_cycle implicit none integer :: i, iostat_val double precision :: T1, T3, rp, gamma_g, cp_g, eta_comp, eta_turb_g double precision :: P_boiler, T_steam, P_cond, eta_turb_s, eta_pump double precision :: T_pinch_min, m_gas, m_steam double precision :: T2s, T2a, T4s, T4a, T_exhaust double precision :: W_comp, W_turb_g, Wnet_brayton, Qin_brayton, eta_brayton double precision :: h1_s, h2s_s, h2a_s, h3_s, h4_s, s1_s double precision :: Wt_s, Wp_s, Wnet_rankine, Qin_rankine, eta_rankine double precision :: Q_HRSG, T_pinch_actual, T_stack double precision :: Wnet_total, Qin_total, eta_combined, heat_rate double precision :: T_sat_boil, hf_boil, hfg_boil, sf_boil, sfg_boil, vf_boil double precision :: T_sat_cond, hf_cond, hfg_cond, sf_cond, sfg_cond, vf_cond double precision :: hg_boil, sg_boil, cp_steam_sh, dT_sh double precision :: rp_i, eta_b_i, eta_c_i, T4a_i double precision :: x4s, h4s_s_val, h4a_s_val double precision, parameter :: Ru = 8.314462d0 read(*,*,iostat=iostat_val) T1 if (iostat_val /= 0) then write(*,*) 'ERROR: Invalid compressor inlet temperature.' stop end if read(*,*,iostat=iostat_val) T3 read(*,*,iostat=iostat_val) rp read(*,*,iostat=iostat_val) gamma_g read(*,*,iostat=iostat_val) cp_g read(*,*,iostat=iostat_val) eta_comp read(*,*,iostat=iostat_val) eta_turb_g read(*,*,iostat=iostat_val) P_boiler read(*,*,iostat=iostat_val) T_steam read(*,*,iostat=iostat_val) P_cond read(*,*,iostat=iostat_val) eta_turb_s read(*,*,iostat=iostat_val) eta_pump read(*,*,iostat=iostat_val) T_pinch_min read(*,*,iostat=iostat_val) m_gas if (iostat_val /= 0) then write(*,*) 'ERROR: Failed to read all combined cycle inputs.' stop end if if (T1<=0.0d0.or.T3<=0.0d0.or.rp<=1.0d0) then write(*,*) 'ERROR: T1, T3 positive and rp > 1 required.' stop end if if (gamma_g<=1.0d0) gamma_g = 1.4d0 if (cp_g<=0.0d0) cp_g = 1005.0d0 if (eta_comp<=0.0d0.or.eta_comp>1.0d0) eta_comp = 0.86d0 if (eta_turb_g<=0.0d0.or.eta_turb_g>1.0d0) eta_turb_g = 0.90d0 if (eta_turb_s<=0.0d0.or.eta_turb_s>1.0d0) eta_turb_s = 0.88d0 if (eta_pump<=0.0d0.or.eta_pump>1.0d0) eta_pump = 0.85d0 if (T_pinch_min<=0.0d0) T_pinch_min = 15.0d0 if (m_gas<=0.0d0) m_gas = 100.0d0 if (P_boiler<=0.0d0) P_boiler = 8.0d6 if (P_cond<=0.0d0) P_cond = 10.0d3 ! ===================== BRAYTON (GAS TURBINE) ===================== ! 1→2: Isentropic compression T2s = T1 * rp**((gamma_g-1.0d0)/gamma_g) T2a = T1 + (T2s - T1) / eta_comp ! 3→4: Isentropic expansion T4s = T3 / rp**((gamma_g-1.0d0)/gamma_g) T4a = T3 - eta_turb_g * (T3 - T4s) W_comp = cp_g * (T2a - T1) ! J/kg gas W_turb_g = cp_g * (T3 - T4a) ! J/kg gas Wnet_brayton = W_turb_g - W_comp Qin_brayton = cp_g * (T3 - T2a) eta_brayton = Wnet_brayton / max(Qin_brayton, 1.0d-10) T_exhaust = T4a ! gas turbine exhaust temperature ! ===================== HRSG ===================== ! Saturation properties at steam boiler and condenser pressures call sat_props(P_boiler, T_sat_boil, hf_boil, hfg_boil, sf_boil, sfg_boil, vf_boil) call sat_props(P_cond, T_sat_cond, hf_cond, hfg_cond, sf_cond, sfg_cond, vf_cond) ! Superheated steam state at turbine inlet hg_boil = hf_boil + hfg_boil sg_boil = sf_boil + sfg_boil cp_steam_sh = 2100.0d0 ! J/(kg K) approximate dT_sh = T_steam - T_sat_boil if (dT_sh < 0.0d0) dT_sh = 0.0d0 h1_s = hg_boil + cp_steam_sh * dT_sh s1_s = sg_boil + cp_steam_sh * log(max(T_steam/T_sat_boil, 1.0d0)) ! Feedwater state (saturated liquid at condenser pressure, pumped) h4_s = hf_cond h2a_s = hf_cond + vf_cond * (P_boiler - P_cond) / eta_pump ! HRSG energy balance: m_gas * cp_g * (T_exhaust - T_stack) = m_steam * (h1_s - h2a_s) ! Pinch point constraint: T_stack >= T_sat_boil + T_pinch_min (simplified) T_pinch_actual = T_exhaust - T_sat_boil T_stack = T_sat_boil + T_pinch_min if (T_stack > T_exhaust - 10.0d0) then T_stack = T_exhaust - 10.0d0 end if if (T_stack < T1 + 20.0d0) T_stack = T1 + 20.0d0 Q_HRSG = m_gas * cp_g * (T_exhaust - T_stack) Qin_rankine = h1_s - h2a_s m_steam = Q_HRSG / max(Qin_rankine, 1.0d-10) ! ===================== RANKINE (STEAM TURBINE) ===================== ! Isentropic expansion to P_cond x4s = (s1_s - sf_cond) / max(sfg_cond, 1.0d-10) if (x4s > 1.0d0) then h4s_s_val = hf_cond + hfg_cond + cp_steam_sh * & (T_sat_cond*exp((s1_s-sg_boil)/cp_steam_sh) - T_sat_cond) else h4s_s_val = hf_cond + x4s * hfg_cond end if h4a_s_val = h1_s - eta_turb_s * (h1_s - h4s_s_val) Wt_s = h1_s - h4a_s_val ! J/kg steam Wp_s = h2a_s - hf_cond ! J/kg steam Wnet_rankine = Wt_s - Wp_s eta_rankine = Wnet_rankine / max(Qin_rankine, 1.0d-10) ! ===================== COMBINED ===================== Wnet_total = m_gas * Wnet_brayton + m_steam * Wnet_rankine ! W total Qin_total = m_gas * Qin_brayton eta_combined = Wnet_total / max(Qin_total, 1.0d-10) heat_rate = 3600.0d0 / max(eta_combined, 1.0d-10) ! kJ/kWh write(*,'(A)') '============================================================' write(*,'(A)') ' COMBINED CYCLE (BRAYTON + RANKINE) ENGINE' write(*,'(A)') '============================================================' write(*,*) write(*,'(A)') '--- BRAYTON (GAS TURBINE) -----------------------------------' write(*,'(A,F12.2,A)') ' T1 (compressor inlet) = ', T1, ' K' write(*,'(A,F12.2,A)') ' T2s (isentropic) = ', T2s, ' K' write(*,'(A,F12.2,A)') ' T2a (actual) = ', T2a, ' K' write(*,'(A,F12.2,A)') ' T3 (turbine inlet) = ', T3, ' K' write(*,'(A,F12.2,A)') ' T4s (isentropic) = ', T4s, ' K' write(*,'(A,F12.2,A)') ' T4a (exhaust) = ', T4a, ' K' write(*,'(A,F12.2)') ' Pressure Ratio rp = ', rp write(*,'(A,F12.2,A)') ' Compressor Work = ', W_comp/1000, ' kJ/kg' write(*,'(A,F12.2,A)') ' Gas Turbine Work = ', W_turb_g/1000, ' kJ/kg' write(*,'(A,F12.2,A)') ' Brayton Net Work = ', Wnet_brayton/1000, ' kJ/kg' write(*,'(A,F12.2,A)') ' Brayton Heat Input = ', Qin_brayton/1000, ' kJ/kg' write(*,'(A,F10.4)') ' Brayton Efficiency = ', eta_brayton write(*,*) write(*,'(A)') '--- HRSG (HEAT RECOVERY) ------------------------------------' write(*,'(A,F12.2,A)') ' Gas Exhaust Temperature = ', T_exhaust, ' K' write(*,'(A,F12.2,A)') ' Stack Temperature = ', T_stack, ' K' write(*,'(A,F12.2,A)') ' Pinch Temperature Delta = ', T_pinch_actual, ' K' write(*,'(A,F12.2,A)') ' HRSG Heat Recovery = ', Q_HRSG/1.0d6, ' MW' write(*,'(A,ES12.4,A)') ' Steam Mass Flow = ', m_steam, ' kg/s' write(*,'(A,F12.4)') ' Steam-to-Gas Mass Ratio = ', m_steam/m_gas write(*,*) write(*,'(A)') '--- RANKINE (STEAM TURBINE) ---------------------------------' write(*,'(A,ES12.4,A)') ' Boiler Pressure = ', P_boiler, ' Pa' write(*,'(A,F12.2,A)') ' Steam Temperature = ', T_steam, ' K' write(*,'(A,ES12.4,A)') ' Condenser Pressure = ', P_cond, ' Pa' write(*,'(A,F12.2,A)') ' Steam Turbine Work = ', Wt_s/1000, ' kJ/kg' write(*,'(A,F12.2,A)') ' Pump Work = ', Wp_s/1000, ' kJ/kg' write(*,'(A,F12.2,A)') ' Rankine Net Work = ', Wnet_rankine/1000, ' kJ/kg' write(*,'(A,F10.4)') ' Rankine Efficiency = ', eta_rankine write(*,*) write(*,'(A)') '--- COMBINED CYCLE PERFORMANCE ------------------------------' write(*,'(A,F12.2,A)') ' Brayton Power = ', m_gas*Wnet_brayton/1.0d6, ' MW' write(*,'(A,F12.2,A)') ' Rankine Power = ', m_steam*Wnet_rankine/1.0d6, ' MW' write(*,'(A,F12.2,A)') ' Total Power = ', Wnet_total/1.0d6, ' MW' write(*,'(A,F12.2,A)') ' Total Heat Input = ', Qin_total/1.0d6, ' MW' write(*,'(A,F10.4)') ' Combined Efficiency = ', eta_combined write(*,'(A,F10.2,A)') ' Combined Efficiency = ', eta_combined*100, ' percent' write(*,'(A,F10.2,A)') ' Heat Rate = ', heat_rate, ' kJ/kWh' write(*,'(A,F10.4)') ' eta_check (1-(1-eB)(1-eR))= ', & 1.0d0-(1.0d0-eta_brayton)*(1.0d0-eta_rankine*Q_HRSG/(m_gas*Qin_brayton)) write(*,*) ! Efficiency vs pressure ratio sweep write(*,'(A)') '--- EFFICIENCY VS PRESSURE RATIO SWEEP ----------------------' write(*,'(A)') ' rp eta_Brayton eta_Combined T_exhaust[K]' write(*,'(A)') ' -----------------------------------------------------------' do i = 1, 50 rp_i = 2.0d0 + 38.0d0*dble(i-1)/49.0d0 T2s = T1 * rp_i**((gamma_g-1.0d0)/gamma_g) T2a = T1 + (T2s - T1)/eta_comp T4s = T3 / rp_i**((gamma_g-1.0d0)/gamma_g) T4a_i = T3 - eta_turb_g*(T3 - T4s) eta_b_i = (cp_g*(T3-T4a_i) - cp_g*(T2a-T1)) / max(cp_g*(T3-T2a), 1.0d-10) ! Simplified combined if (T4a_i > T_stack) then eta_c_i = eta_b_i + (1.0d0-eta_b_i)*eta_rankine* & (T4a_i-T_stack)/(T4a_i-T1) else eta_c_i = eta_b_i end if write(*,'(F10.2,2X,F10.5,2X,F10.5,2X,F12.2)') rp_i, eta_b_i, eta_c_i, T4a_i end do write(*,*) ! T-Q diagram for HRSG write(*,'(A)') '--- HRSG TQ DIAGRAM DATA ------------------------------------' write(*,'(A)') ' Q_frac T_gas[K] T_steam[K]' write(*,'(A)') ' -------------------------------------------' do i = 0, 50 call hrsg_tq(dble(i)/50.0d0, T_exhaust, T_stack, & T_sat_boil, T_steam, hf_cond, h2a_s, h1_s, hf_boil, hg_boil) end do write(*,*) write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------' write(*,'(A)') ' Brayton: T2s=T1*rp^((g-1)/g); T4s=T3/rp^((g-1)/g).' write(*,'(A)') ' eta_combined ~ 1-(1-eta_B)(1-eta_R * Q_HRSG/Qin_B).' write(*,'(A)') ' HRSG: m_gas*cp*(T_exh-T_stack)=m_steam*(h1-h2a).' 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) Tsat = 373.15d0 + 42.0d0*logP - 0.8d0*logP**2 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 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 if (sf < 0.3d3) sf = 0.3d3 + 0.5d3*P_MPa sfg = hfg / Tsat vf = 0.001d0 * (1.0d0 + 0.0002d0*P_MPa) end subroutine sat_props subroutine hrsg_tq(frac, T_exh, T_stk, T_sat, T_sh, hfc, h2a, h1, hfb, hgb) implicit none double precision, intent(in) :: frac, T_exh, T_stk, T_sat, T_sh double precision, intent(in) :: hfc, h2a, h1, hfb, hgb double precision :: Q_tot, Q_i, T_gas, T_stm double precision :: Q_eco, Q_evap, Q_sh, frac_eco, frac_evap Q_tot = h1 - h2a Q_eco = hfb - h2a Q_evap = hgb - hfb Q_sh = h1 - hgb frac_eco = Q_eco / max(Q_tot, 1.0d-10) frac_evap = Q_evap / max(Q_tot, 1.0d-10) Q_i = frac * Q_tot T_gas = T_exh - (T_exh - T_stk)*frac if (frac < frac_eco) then T_stm = h2a/max(hfb,1.0d-10)*T_sat * frac/max(frac_eco,1.0d-10) T_stm = T_sat * frac/max(frac_eco,1.0d-10) + 273.15d0*(1.0d0-frac/max(frac_eco,1.0d-10)) if (T_stm < 300.0d0) T_stm = 300.0d0 else if (frac < frac_eco + frac_evap) then T_stm = T_sat else T_stm = T_sat + (T_sh - T_sat)*(frac - frac_eco - frac_evap)/max(1.0d0-frac_eco-frac_evap,1.0d-10) end if write(*,'(F10.4,2X,F12.2,2X,F12.2)') frac, T_gas, T_stm end subroutine hrsg_tq end program combined_cycle