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