💻 Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Combined Cycle Power Plant
Core Numerical Engine in Fortran 90 • 29 total downloads
combined_cycle.f90
! =========================================================================
! Source File: combined_cycle.f90
! =========================================================================
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
Solver Description
Model gas-steam combined cycle power plants. Integrates a gas turbine Brayton topping cycle with a steam Rankine bottoming cycle. Evaluates turbine exhaust heat recovery, performs pinch point temperature analysis in the Heat Recovery Steam Generator (HRSG) to determine steam mass flow rates, and calculates Brayton, Rankine, and overall combined plant thermal efficiencies and power splits.
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
🛠️ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -O3 combined_cycle.f90 -o combined_cycle
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
combined_cycle < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Compressor inlet temperature T1 [K]
288.15
! Turbine inlet temperature T3 [K]
1573.15
! Pressure ratio rp
18.0
! Specific heat ratio gamma (gas)
1.4
! Specific heat cp [J/kg-K] (gas)
1005.0
! Compressor isentropic efficiency
0.86
! Gas turbine isentropic efficiency
0.90
! Steam boiler pressure Pb [Pa]
8.0e6
! Steam temperature Ts [K]
813.15
! Steam condenser pressure Pc [Pa]
10.0e3
! Steam turbine isentropic efficiency
0.88
! Pump isentropic efficiency
0.85
! Min pinch temp delta [K]
15.0
! Gas mass flow rate [kg/s]
400.0
288.15
! Turbine inlet temperature T3 [K]
1573.15
! Pressure ratio rp
18.0
! Specific heat ratio gamma (gas)
1.4
! Specific heat cp [J/kg-K] (gas)
1005.0
! Compressor isentropic efficiency
0.86
! Gas turbine isentropic efficiency
0.90
! Steam boiler pressure Pb [Pa]
8.0e6
! Steam temperature Ts [K]
813.15
! Steam condenser pressure Pc [Pa]
10.0e3
! Steam turbine isentropic efficiency
0.88
! Pump isentropic efficiency
0.85
! Min pinch temp delta [K]
15.0
! Gas mass flow rate [kg/s]
400.0