🔬
Solver Purpose & Physical Scope
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.
📂 Discipline: Thermo
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 284 times
📄 Source File:
rankine_orc.f90
📁 calcul/Thermodynamics /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 rankine_orc.f90 -o rankine_orc
2. Execution with input.txt redirection:
rankine_orc < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 423.15 303.15 5.0 0.82 0.75 15.0 10.0
Parameter Description:
Organic fluid (1=R245fa, 2=R134a, 3=n-Pentane, 4=Toluene, 5=R1233zd) Heat source temperature Tsource [K] Cooling sink temperature Tsink [K] Superheat delta Tsh [K] Turbine isentropic efficiency Pump isentropic efficiency Mass flow rate [kg/s] Evaporator pinch point delta T [K]