program brayton_cycle implicit none ! Constants for air (cold-air-standard assumptions) double precision, parameter :: Cp = 1.005d0 ! kJ/(kg.K) double precision, parameter :: Cv = 0.718d0 ! kJ/(kg.K) double precision, parameter :: R = 0.287d0 ! kJ/(kg.K) double precision, parameter :: gamma = 1.4d0 double precision, parameter :: T_ref = 273.15d0 ! K for entropy double precision, parameter :: P_ref = 100.0d0 ! kPa for entropy ! Inputs double precision :: T1_c, P1, r_p, T3_c, eta_c, eta_t, m_dot integer :: opt_regen, opt_inter, opt_reheat double precision :: eta_regen ! State properties (Temperatures in K, pressures in kPa, entropies in kJ/kgK) double precision :: T1, P1_val, s1 double precision :: T1a, P1a, s1a, T1as double precision :: T1b, P1b, s1b double precision :: T2, P2, s2, T2s double precision :: T5, P5, s5 double precision :: T3, P3, s3 double precision :: T3a, P3a, s3a, T3as double precision :: T3b, P3b, s3b double precision :: T4, P4, s4, T4s double precision :: T6, P6, s6 ! Cycle metrics double precision :: W_c1, W_c2, W_c, W_t1, W_t2, W_t, W_net double precision :: Q_in, Q_out, eta_thermal, BWR double precision :: W_net_total, Q_in_total ! Read from stdin read(*,*) T1_c read(*,*) P1 read(*,*) r_p read(*,*) T3_c read(*,*) eta_c read(*,*) eta_t read(*,*) m_dot read(*,*) opt_regen read(*,*) eta_regen read(*,*) opt_inter read(*,*) opt_reheat ! Convert temperatures to Kelvin T1 = T1_c + 273.15d0 T3 = T3_c + 273.15d0 P1_val = P1 ! ------------------------------------------------------------- ! COMPRESSION PROCESS ! ------------------------------------------------------------- s1 = entropy(T1, P1_val) if (opt_inter == 1) then ! 2-Stage Compression with Perfect Intercooling (intermediate ratio = sqrt(r_p)) P1a = P1_val * sqrt(r_p) T1as = T1 * (P1a / P1_val)**((gamma - 1.0d0) / gamma) T1a = T1 + (T1as - T1) / eta_c s1a = entropy(T1a, P1a) W_c1 = Cp * (T1a - T1) ! Intercooling back to T1 T1b = T1 P1b = P1a s1b = entropy(T1b, P1b) ! Second compression stage P2 = P1_val * r_p T2s = T1b * (P2 / P1b)**((gamma - 1.0d0) / gamma) T2 = T1b + (T2s - T1b) / eta_c s2 = entropy(T2, P2) W_c2 = Cp * (T2 - T1b) W_c = W_c1 + W_c2 else ! Single Stage Compression P2 = P1_val * r_p T2s = T1 * (P2 / P1_val)**((gamma - 1.0d0) / gamma) T2 = T1 + (T2s - T1) / eta_c s2 = entropy(T2, P2) W_c = Cp * (T2 - T1) ! Fill intermediate variables to avoid unitialized references T1a = T2 P1a = P2 s1a = s2 T1b = T2 P1b = P2 s1b = s2 end if ! ------------------------------------------------------------- ! TURBINE EXPANSION PROCESS ! ------------------------------------------------------------- P3 = P2 s3 = entropy(T3, P3) if (opt_reheat == 1) then ! 2-Stage Expansion with Perfect Reheat P3a = P1_val * sqrt(r_p) T3as = T3 / (P3 / P3a)**((gamma - 1.0d0) / gamma) T3a = T3 - eta_t * (T3 - T3as) s3a = entropy(T3a, P3a) W_t1 = Cp * (T3 - T3a) ! Reheat back to T3 T3b = T3 P3b = P3a s3b = entropy(T3b, P3b) ! Second expansion stage P4 = P1_val T4s = T3b / (P3b / P4)**((gamma - 1.0d0) / gamma) T4 = T3b - eta_t * (T3b - T4s) s4 = entropy(T4, P4) W_t2 = Cp * (T3b - T4) W_t = W_t1 + W_t2 else ! Single Stage Expansion P4 = P1_val T4s = T3 / (P3 / P4)**((gamma - 1.0d0) / gamma) T4 = T3 - eta_t * (T3 - T4s) s4 = entropy(T4, P4) W_t = Cp * (T3 - T4) ! Fill intermediate variables T3a = T4 P3a = P4 s3a = s4 T3b = T4 P3b = P4 s3b = s4 end if ! ------------------------------------------------------------- ! REGENERATION ! ------------------------------------------------------------- P5 = P2 P6 = P4 if (opt_regen == 1 .and. T4 > T2) then T5 = T2 + eta_regen * (T4 - T2) T6 = T4 - eta_regen * (T4 - T2) else T5 = T2 T6 = T4 end if s5 = entropy(T5, P5) s6 = entropy(T6, P6) ! ------------------------------------------------------------- ! PERFORMANCE METRICS ! ------------------------------------------------------------- W_net = W_t - W_c ! Heat input: Combustor (5 -> 3) + optional Reheater (3a -> 3b) Q_in = Cp * (T3 - T5) if (opt_reheat == 1) then Q_in = Q_in + Cp * (T3b - T3a) end if ! Heat rejected: Exhaust (6 -> 1) + optional Intercooler (1a -> 1b) Q_out = Cp * (T6 - T1) if (opt_inter == 1) then Q_out = Q_out + Cp * (T1a - T1b) end if eta_thermal = W_net / Q_in BWR = W_c / W_t W_net_total = m_dot * W_net Q_in_total = m_dot * Q_in ! ------------------------------------------------------------- ! OUTPUT IN PARSABLE FORMAT ! ------------------------------------------------------------- write(*,'(A,F14.4)') 'T1=', T1_c write(*,'(A,F14.4)') 'P1=', P1 write(*,'(A,F14.4)') 'T2=', T2 - 273.15d0 write(*,'(A,F14.4)') 'P2=', P2 write(*,'(A,F14.4)') 'T3=', T3_c write(*,'(A,F14.4)') 'P3=', P3 write(*,'(A,F14.4)') 'T4=', T4 - 273.15d0 write(*,'(A,F14.4)') 'P4=', P4 write(*,'(A,F14.4)') 'T5=', T5 - 273.15d0 write(*,'(A,F14.4)') 'T6=', T6 - 273.15d0 write(*,'(A,F14.4)') 'T1A=', T1a - 273.15d0 write(*,'(A,F14.4)') 'P1A=', P1a write(*,'(A,F14.4)') 'T1B=', T1b - 273.15d0 write(*,'(A,F14.4)') 'P1B=', P1b write(*,'(A,F14.4)') 'T3A=', T3a - 273.15d0 write(*,'(A,F14.4)') 'P3A=', P3a write(*,'(A,F14.4)') 'T3B=', T3b - 273.15d0 write(*,'(A,F14.4)') 'P3B=', P3b write(*,'(A,F14.6)') 'S1=', s1 write(*,'(A,F14.6)') 'S2=', s2 write(*,'(A,F14.6)') 'S3=', s3 write(*,'(A,F14.6)') 'S4=', s4 write(*,'(A,F14.6)') 'S5=', s5 write(*,'(A,F14.6)') 'S6=', s6 write(*,'(A,F14.4)') 'WC=', W_c write(*,'(A,F14.4)') 'WT=', W_t write(*,'(A,F14.4)') 'WNET=', W_net write(*,'(A,F14.4)') 'QIN=', Q_in write(*,'(A,F14.4)') 'QOUT=', Q_out write(*,'(A,F14.4)') 'ETA_TH=', eta_thermal * 100.0d0 write(*,'(A,F14.4)') 'BWR=', BWR * 100.0d0 write(*,'(A,F14.2)') 'WNET_TOT=', W_net_total write(*,'(A,F14.2)') 'QIN_TOT=', Q_in_total contains ! Simple entropy calculation function double precision function entropy(T_k, P_kpa) double precision, intent(in) :: T_k, P_kpa entropy = Cp * log(T_k / T_ref) - R * log(P_kpa / P_ref) end function entropy end program brayton_cycle