💻 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.
Organic Rankine Cycle (ORC)
Core Numerical Engine in Fortran 90 • 27 total downloads
rankine_orc.f90
! =========================================================================
! Source File: rankine_orc.f90
! =========================================================================
program rankine_orc
implicit none
integer :: fluid_type, iostat_val, i, n_sweep, iter
double precision :: T_source, T_sink, T_superheat, eta_t, eta_p, mdot, T_pinch
double precision :: Tc, Pc, Tb, M_fluid, cp_liq, cp_vap, hfg_Tb, rho_liq
double precision :: T_evap, T_cond, P_evap, P_cond
double precision :: h1, s1, h2s, h2a, h3, h4, v_f
double precision :: T2s, hfg_evap, hfg_cond
double precision :: Wt, Wp, Wnet, Qin, Qout, eta_th, eta_Carnot, eta_rel
double precision :: P_output, SSC, bwr
double precision :: Ts_sw, eta_sw, Wnet_sw
double precision :: Te_sw, Tc_sw, Pe_sw, Pcd_sw
double precision :: h1s, s1s, h2ss, h2as, h3s, h4s, vfs
double precision :: hfg_e_sw, hfg_c_sw, T2ss, Qins
character(len=40) :: fluid_name
read(*,*,iostat=iostat_val) fluid_type
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid fluid type.'; stop; end if
read(*,*,iostat=iostat_val) T_source
read(*,*,iostat=iostat_val) T_sink
read(*,*,iostat=iostat_val) T_superheat
read(*,*,iostat=iostat_val) eta_t
read(*,*,iostat=iostat_val) eta_p
read(*,*,iostat=iostat_val) mdot
read(*,*,iostat=iostat_val) T_pinch
if (iostat_val /= 0) then; write(*,*) 'ERROR: Failed to read all inputs.'; stop; end if
if (T_source <= 0.0d0) T_source = 423.15d0
if (T_sink <= 0.0d0) T_sink = 303.15d0
if (T_superheat < 0.0d0) T_superheat = 5.0d0
if (eta_t <= 0.0d0 .or. eta_t > 1.0d0) eta_t = 0.82d0
if (eta_p <= 0.0d0 .or. eta_p > 1.0d0) eta_p = 0.75d0
if (mdot <= 0.0d0) mdot = 10.0d0
if (T_pinch <= 0.0d0) T_pinch = 10.0d0
if (T_source <= T_sink + 2.0d0*T_pinch) then
write(*,*) 'ERROR: T_source must exceed T_sink + 2*T_pinch.'
stop
end if
! ── Fluid properties ───────────────────────────────────────
select case(fluid_type)
case(1)
fluid_name='R245fa (1,1,1,3,3-Pentafluoropropane)'
Tc=427.2d0; Pc=3.651d0; Tb=288.3d0; M_fluid=134.05d0
cp_liq=1360.0d0; cp_vap=900.0d0; hfg_Tb=196000.0d0; rho_liq=1320.0d0
case(2)
fluid_name='R134a (1,1,1,2-Tetrafluoroethane)'
Tc=374.2d0; Pc=4.059d0; Tb=247.1d0; M_fluid=102.03d0
cp_liq=1430.0d0; cp_vap=850.0d0; hfg_Tb=217000.0d0; rho_liq=1210.0d0
case(3)
fluid_name='n-Pentane (C5H12)'
Tc=469.7d0; Pc=3.370d0; Tb=309.2d0; M_fluid=72.15d0
cp_liq=2320.0d0; cp_vap=1670.0d0; hfg_Tb=358000.0d0; rho_liq=626.0d0
case(4)
fluid_name='Toluene (C7H8)'
Tc=591.8d0; Pc=4.109d0; Tb=383.8d0; M_fluid=92.14d0
cp_liq=1710.0d0; cp_vap=1130.0d0; hfg_Tb=363000.0d0; rho_liq=867.0d0
case default
fluid_name='R1233zd(E) (trans-1-Chloro-3,3,3-trifluoropropene)'
Tc=438.8d0; Pc=3.624d0; Tb=291.4d0; M_fluid=130.50d0
cp_liq=1250.0d0; cp_vap=820.0d0; hfg_Tb=195000.0d0; rho_liq=1290.0d0
fluid_type=5
end select
! ── Evaporator / condenser conditions ──────────────────────
T_evap = T_source - T_pinch
T_cond = T_sink + T_pinch
if (T_evap >= Tc) T_evap = Tc - 5.0d0
if (T_cond >= T_evap) then
write(*,*) 'ERROR: T_evap must exceed T_cond.'
stop
end if
! Pressure from Clausius-Clapeyron approximation
! ln(P/Pc) ~ (hfg*M/(R*1000)) * (1/Tc - 1/T) simplified
! Use: P = Pc * exp( -5.0*(Tb/T - 1.0)*(Tc/Tb) ) (rough fit)
P_evap = Pc * exp(-5.0d0 * (Tb/T_evap - 1.0d0) * Tc/Tb)
if (P_evap > 0.90d0*Pc) P_evap = 0.90d0*Pc
if (P_evap < 0.01d0) P_evap = 0.01d0
P_cond = Pc * exp(-5.0d0 * (Tb/T_cond - 1.0d0) * Tc/Tb)
if (P_cond > P_evap) P_cond = P_evap * 0.1d0
if (P_cond < 0.001d0) P_cond = 0.001d0
! Latent heat adjusted with temperature (Trouton-Watson)
hfg_evap = hfg_Tb * ((Tc - T_evap)/(Tc - Tb))**0.38d0
hfg_cond = hfg_Tb * ((Tc - T_cond)/(Tc - Tb))**0.38d0
if (hfg_evap < 0.0d0) hfg_evap = hfg_Tb*0.5d0
if (hfg_cond < 0.0d0) hfg_cond = hfg_Tb*0.5d0
! ── State points (J/kg, J/(kg.K)) ─────────────────────────
! Reference: saturated liquid at Tb => h=0, s=0
! State 1: superheated vapor at P_evap
h1 = cp_liq*(T_evap - Tb) + hfg_evap + cp_vap*T_superheat
s1 = cp_liq*log(T_evap/Tb) + hfg_evap/T_evap + cp_vap*log((T_evap+T_superheat)/T_evap)
! State 2s: isentropic expansion to P_cond
! For dry ORC fluid, exit is superheated
! s2s = s1 => cp_liq*ln(T_cond/Tb) + hfg_cond/T_cond + cp_vap*ln(T2s/T_cond) = s1
! => T2s = T_cond * exp( (s1 - cp_liq*ln(T_cond/Tb) - hfg_cond/T_cond) / cp_vap )
T2s = T_cond * exp((s1 - cp_liq*log(T_cond/Tb) - hfg_cond/T_cond) / cp_vap)
if (T2s < T_cond) T2s = T_cond
h2s = cp_liq*(T_cond - Tb) + hfg_cond + cp_vap*(T2s - T_cond)
! State 2a: actual turbine exit
h2a = h1 - eta_t * (h1 - h2s)
! State 3: saturated liquid at condenser
h3 = cp_liq * (T_cond - Tb)
! State 4: compressed liquid after pump
v_f = 1.0d0 / rho_liq ! m3/kg
h4 = h3 + v_f * (P_evap - P_cond) * 1.0d6 / eta_p ! Pa*m3/kg = J/kg
! ── Energy analysis (kJ/kg) ────────────────────────────────
Wt = (h1 - h2a) / 1000.0d0
Wp = (h4 - h3) / 1000.0d0
Wnet = Wt - Wp
Qin = (h1 - h4) / 1000.0d0
Qout = (h2a - h3) / 1000.0d0
eta_th = Wnet / max(Qin, 1.0d-10)
eta_Carnot = 1.0d0 - T_sink / T_source
eta_rel = eta_th / max(eta_Carnot, 1.0d-10)
P_output = mdot * Wnet
bwr = Wp / max(Wt, 1.0d-10)
SSC = 3600.0d0 / max(Wnet, 1.0d-10)
! ── Output ──────────────────────────────────────────────────
write(*,'(A)') '============================================================'
write(*,'(A)') ' ORGANIC RANKINE CYCLE (ORC)'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUTS --------------------------------------------------'
write(*,'(A,A)') ' Working Fluid = ', trim(fluid_name)
write(*,'(A,F12.2,A)') ' Heat Source Temperature = ', T_source, ' K'
write(*,'(A,F12.2,A)') ' Heat Sink Temperature = ', T_sink, ' K'
write(*,'(A,F12.2,A)') ' Superheat = ', T_superheat, ' K'
write(*,'(A,F10.4)') ' Turbine Isentropic Eff = ', eta_t
write(*,'(A,F10.4)') ' Pump Isentropic Eff = ', eta_p
write(*,'(A,F12.4,A)') ' Mass Flow Rate = ', mdot, ' kg/s'
write(*,'(A,F12.2,A)') ' Pinch Point DeltaT = ', T_pinch, ' K'
write(*,*)
write(*,'(A)') '--- FLUID PROPERTIES ----------------------------------------'
write(*,'(A,F12.2,A)') ' Critical Temperature Tc = ', Tc, ' K'
write(*,'(A,F12.4,A)') ' Critical Pressure Pc = ', Pc, ' MPa'
write(*,'(A,F12.2,A)') ' Normal Boiling Point Tb = ', Tb, ' K'
write(*,'(A,F12.2,A)') ' Molar Mass M = ', M_fluid, ' g/mol'
write(*,'(A,F12.2,A)') ' cp_liquid = ', cp_liq, ' J/(kg.K)'
write(*,'(A,F12.2,A)') ' cp_vapor = ', cp_vap, ' J/(kg.K)'
write(*,'(A,F12.2,A)') ' hfg at Tb = ', hfg_Tb/1000.0d0, ' kJ/kg'
write(*,'(A,F12.2,A)') ' Liquid density = ', rho_liq, ' kg/m3'
write(*,*)
write(*,'(A)') '--- CYCLE PRESSURES & TEMPERATURES --------------------------'
write(*,'(A,F12.2,A)') ' T_evaporator = ', T_evap, ' K'
write(*,'(A,F12.2,A)') ' T_condenser = ', T_cond, ' K'
write(*,'(A,F12.6,A)') ' P_evaporator = ', P_evap, ' MPa'
write(*,'(A,F12.6,A)') ' P_condenser = ', P_cond, ' MPa'
write(*,'(A,F12.2,A)') ' hfg at T_evap = ', hfg_evap/1000.0d0, ' kJ/kg'
write(*,'(A,F12.2,A)') ' hfg at T_cond = ', hfg_cond/1000.0d0, ' kJ/kg'
write(*,*)
write(*,'(A)') '--- STATE POINTS (kJ/kg, relative to sat liq at Tb) --------'
write(*,'(A,F12.2,A)') ' h1 (turb inlet) = ', h1/1000.0d0, ' kJ/kg'
write(*,'(A,F12.2,A)') ' h2s (turb exit, isen) = ', h2s/1000.0d0, ' kJ/kg'
write(*,'(A,F12.2,A)') ' h2a (turb exit, actual) = ', h2a/1000.0d0, ' kJ/kg'
write(*,'(A,F12.2,A)') ' h3 (cond exit, sat liq) = ', h3/1000.0d0, ' kJ/kg'
write(*,'(A,F12.2,A)') ' h4 (pump exit) = ', h4/1000.0d0, ' kJ/kg'
write(*,'(A,F12.2,A)') ' T2s (turb exit temp) = ', T2s, ' K'
write(*,'(A,F12.2,A)') ' T1 (turb inlet temp) = ', T_evap+T_superheat, ' K'
write(*,*)
write(*,'(A)') '--- ENERGY PER UNIT MASS ------------------------------------'
write(*,'(A,F12.4,A)') ' Turbine Work Wt = ', Wt, ' kJ/kg'
write(*,'(A,F12.4,A)') ' Pump Work Wp = ', Wp, ' kJ/kg'
write(*,'(A,F12.4,A)') ' Net Work Wnet = ', Wnet, ' kJ/kg'
write(*,'(A,F12.4,A)') ' Heat Input Qin = ', Qin, ' kJ/kg'
write(*,'(A,F12.4,A)') ' Heat Rejected Qout = ', Qout, ' kJ/kg'
write(*,*)
write(*,'(A)') '--- CYCLE PERFORMANCE ---------------------------------------'
write(*,'(A,F10.6)') ' Thermal Efficiency = ', eta_th
write(*,'(A,F10.2,A)') ' Thermal Efficiency = ', eta_th*100.0d0, ' percent'
write(*,'(A,F10.6)') ' Carnot Efficiency = ', eta_Carnot
write(*,'(A,F10.2,A)') ' Relative Efficiency = ', eta_rel*100.0d0, ' percent of Carnot'
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 Consump = ', SSC, ' kg/kWh'
write(*,*)
! ── Sensitivity sweep: T_source ────────────────────────────
n_sweep = 40
write(*,'(A)') '--- SENSITIVITY: EFFICIENCY VS SOURCE TEMPERATURE ----------'
write(*,'(A)') ' T_source[K] eta_th Wnet[kJ/kg] P_out[kW]'
write(*,'(A)') ' -----------------------------------------------------------'
do i = 1, n_sweep
Ts_sw = (T_sink+30.0d0) + dble(i-1)*(T_source+80.0d0 - (T_sink+30.0d0))/dble(n_sweep-1)
Te_sw = Ts_sw - T_pinch
if (Te_sw >= Tc) Te_sw = Tc - 5.0d0
Tc_sw = T_cond
Pe_sw = Pc * exp(-5.0d0*(Tb/Te_sw - 1.0d0)*Tc/Tb)
if (Pe_sw > 0.90d0*Pc) Pe_sw = 0.90d0*Pc
if (Pe_sw < 0.01d0) Pe_sw = 0.01d0
Pcd_sw = P_cond
hfg_e_sw = hfg_Tb * ((Tc-Te_sw)/(Tc-Tb))**0.38d0
if (hfg_e_sw < 0.0d0) hfg_e_sw = hfg_Tb*0.3d0
hfg_c_sw = hfg_cond
h1s = cp_liq*(Te_sw-Tb)+hfg_e_sw+cp_vap*T_superheat
s1s = cp_liq*log(Te_sw/Tb)+hfg_e_sw/Te_sw+cp_vap*log((Te_sw+T_superheat)/Te_sw)
T2ss = Tc_sw*exp((s1s-cp_liq*log(Tc_sw/Tb)-hfg_c_sw/Tc_sw)/cp_vap)
if (T2ss < Tc_sw) T2ss = Tc_sw
h2ss = cp_liq*(Tc_sw-Tb)+hfg_c_sw+cp_vap*(T2ss-Tc_sw)
h2as = h1s - eta_t*(h1s-h2ss)
h3s = cp_liq*(Tc_sw-Tb)
vfs = 1.0d0/rho_liq
h4s = h3s + vfs*(Pe_sw-Pcd_sw)*1.0d6/eta_p
Qins = (h1s-h4s)/1000.0d0
Wnet_sw = ((h1s-h2as)-(h4s-h3s))/1000.0d0
if (Qins > 0.0d0) then
eta_sw = Wnet_sw / Qins
else
eta_sw = 0.0d0
end if
write(*,'(F10.2,4X,F10.6,4X,F12.4,4X,F12.2)') Ts_sw, eta_sw, Wnet_sw, mdot*Wnet_sw
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' Pressure: Clausius-Clapeyron approx with empirical fit.'
write(*,'(A)') ' Latent heat: Watson correlation hfg ~ hfg_Tb*((Tc-T)/(Tc-Tb))^0.38.'
write(*,'(A)') ' Isentropic expansion assuming dry fluid (superheated exit).'
write(*,'(A)') ' Pump work: w_p = v_f*(P_evap-P_cond)/eta_p.'
write(*,'(A)') ' Simplified constant cp model (educational accuracy).'
end program rankine_orc
Solver Description
Model Organic Rankine Cycles (ORC) designed for low-temperature waste heat recovery, solar thermal, and geothermal power systems. Computes thermodynamic state points for organic working fluids (R134a, R245fa, n-Pentane, Toluene, R1233zd) using specialized equation-of-state property fits. Evaluates thermal efficiency, net power output, back-work ratio, specific steam consumption (SSC), and pinch-point constraint metrics.
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 rankine_orc.f90 -o rankine_orc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
rankine_orc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Organic fluid (1=R245fa, 2=R134a, 3=n-Pentane, 4=Toluene, 5=R1233zd)
1
! Heat source temperature Tsource [K]
423.15
! Cooling sink temperature Tsink [K]
303.15
! Superheat delta Tsh [K]
5.0
! Turbine isentropic efficiency
0.82
! Pump isentropic efficiency
0.75
! Mass flow rate [kg/s]
15.0
! Evaporator pinch point delta T [K]
10.0
1
! Heat source temperature Tsource [K]
423.15
! Cooling sink temperature Tsink [K]
303.15
! Superheat delta Tsh [K]
5.0
! Turbine isentropic efficiency
0.82
! Pump isentropic efficiency
0.75
! Mass flow rate [kg/s]
15.0
! Evaporator pinch point delta T [K]
10.0