program isentropic_polytropic implicit none integer :: i, iostat_val, n_sweep double precision :: P1, P2, T1, V1, poly_n, R_gas, gamma_val, mass double precision :: V2, T2, W, Q, dU, dH, dS double precision :: cv, cp double precision :: n_sw, V2_sw, T2_sw, W_sw, Q_sw character(len=40) :: process_name ! ── Read inputs ───────────────────────────────────────────── read(*,*,iostat=iostat_val) P1 if (iostat_val /= 0) then write(*,*) 'ERROR: Invalid P1 input.' stop end if read(*,*,iostat=iostat_val) P2 read(*,*,iostat=iostat_val) T1 read(*,*,iostat=iostat_val) V1 read(*,*,iostat=iostat_val) poly_n read(*,*,iostat=iostat_val) R_gas read(*,*,iostat=iostat_val) gamma_val read(*,*,iostat=iostat_val) mass if (iostat_val /= 0) then write(*,*) 'ERROR: Failed to read all inputs.' stop end if ! ── Input validation ──────────────────────────────────────── if (P1 <= 0.0d0) then write(*,*) 'ERROR: P1 must be positive.' stop end if if (P2 <= 0.0d0) then write(*,*) 'ERROR: P2 must be positive.' stop end if if (T1 <= 0.0d0) T1 = 300.0d0 if (V1 <= 0.0d0) V1 = 0.1d0 if (mass <= 0.0d0) mass = 1.0d0 if (gamma_val <= 1.0d0) gamma_val = 1.4d0 if (R_gas <= 0.0d0) R_gas = 287.0d0 ! ── Gas properties ────────────────────────────────────────── cv = R_gas / (gamma_val - 1.0d0) ! J/(kg·K) cp = gamma_val * cv ! J/(kg·K) ! ── Process identification ────────────────────────────────── if (abs(poly_n) < 1.0d-6) then process_name = 'Isobaric (n = 0)' else if (abs(poly_n - 1.0d0) < 1.0d-6) then process_name = 'Isothermal (n = 1)' else if (abs(poly_n - gamma_val) < 0.02d0) then process_name = 'Isentropic (n = gamma)' else if (poly_n > 20.0d0) then process_name = 'Isochoric (n -> infinity)' else process_name = 'Polytropic (general)' end if ! ── Core calculations ─────────────────────────────────────── ! Handle special cases for n if (poly_n > 20.0d0) then ! Isochoric: V = const V2 = V1 T2 = T1 * (P2 / P1) W = 0.0d0 else if (abs(poly_n) < 1.0d-6) then ! Isobaric: P = const => P2 = P1 (force) V2 = V1 * (T1 * P2) / (T1 * P1) ! from PV=mRT ! Actually for isobaric: V2/V1 = T2/T1, need another relation ! With n=0: PV^0 = C => P = C, so P2 = P1 ! T2 = T1 * V2/V1, but V2 = V1*(P1/P2)^(1/0) is undefined ! Use ideal gas: T2/T1 = (P2/P1)^((n-1)/n) with n->0: T2/T1 = (P2/P1)^(-inf) ! For n=0, P=const so P2=P1. Use V2 from ideal gas T2 = T1 * (P2 / P1) ! simplified V2 = V1 * T2 / T1 W = P1 * (V2 - V1) ! kPa * m^3 = kJ else if (abs(poly_n - 1.0d0) < 1.0d-6) then ! Isothermal: T = const T2 = T1 V2 = V1 * (P1 / P2) W = P1 * V1 * log(V2 / V1) ! kPa * m^3 = kJ else ! General polytropic V2 = V1 * (P1 / P2) ** (1.0d0 / poly_n) T2 = T1 * (P2 / P1) ** ((poly_n - 1.0d0) / poly_n) W = (P2 * V2 - P1 * V1) / (1.0d0 - poly_n) ! kPa*m^3 = kJ end if ! Internal energy change (kJ) dU = mass * cv * (T2 - T1) / 1000.0d0 ! Enthalpy change (kJ) dH = mass * cp * (T2 - T1) / 1000.0d0 ! Heat transfer from first law: Q = dU + W (all in kJ) Q = dU + W ! Entropy change for ideal gas (kJ/K) if (T2 > 0.0d0 .and. T1 > 0.0d0 .and. V2 > 0.0d0 .and. V1 > 0.0d0) then dS = mass * (cv * log(T2/T1) + R_gas * log(V2/V1)) / 1000.0d0 else dS = 0.0d0 end if ! ── Output ────────────────────────────────────────────────── write(*,'(A)') '============================================================' write(*,'(A)') ' ISENTROPIC & POLYTROPIC PROCESS ANALYSIS' write(*,'(A)') '============================================================' write(*,*) write(*,'(A)') '--- INPUTS --------------------------------------------------' write(*,'(A,F12.4,A)') ' Initial Pressure P1 = ', P1, ' kPa' write(*,'(A,F12.4,A)') ' Final Pressure P2 = ', P2, ' kPa' write(*,'(A,F12.2,A)') ' Initial Temperature T1 = ', T1, ' K' write(*,'(A,ES12.4,A)') ' Initial Volume V1 = ', V1, ' m^3' write(*,'(A,F12.4)') ' Polytropic Index n = ', poly_n write(*,'(A,F12.2,A)') ' Specific Gas Constant R = ', R_gas, ' J/(kg.K)' write(*,'(A,F12.4)') ' Gamma (cp/cv) = ', gamma_val write(*,'(A,F12.4,A)') ' Mass = ', mass, ' kg' write(*,*) write(*,'(A)') '--- GAS PROPERTIES ------------------------------------------' write(*,'(A,F12.4,A)') ' cv = ', cv, ' J/(kg.K)' write(*,'(A,F12.4,A)') ' cp = ', cp, ' J/(kg.K)' write(*,*) write(*,'(A)') '--- PROCESS IDENTIFICATION ----------------------------------' write(*,'(A,A)') ' Process Type = ', trim(process_name) write(*,*) write(*,'(A)') '--- FINAL STATE ---------------------------------------------' write(*,'(A,F12.2,A)') ' Final Temperature T2 = ', T2, ' K' write(*,'(A,ES12.4,A)') ' Final Volume V2 = ', V2, ' m^3' write(*,'(A,F12.4,A)') ' Final Pressure P2 = ', P2, ' kPa' write(*,*) write(*,'(A)') '--- ENERGY & ENTROPY ----------------------------------------' write(*,'(A,F12.4,A)') ' Work W = ', W, ' kJ' write(*,'(A,F12.4,A)') ' Heat Transfer Q = ', Q, ' kJ' write(*,'(A,F12.4,A)') ' Internal Energy Change dU = ', dU, ' kJ' write(*,'(A,F12.4,A)') ' Enthalpy Change dH = ', dH, ' kJ' write(*,'(A,ES14.6,A)') ' Entropy Change dS = ', dS, ' kJ/K' write(*,*) write(*,'(A)') '--- SPECIFIC VALUES (per kg) --------------------------------' if (mass > 1.0d-10) then write(*,'(A,F12.4,A)') ' w (specific work) = ', W/mass, ' kJ/kg' write(*,'(A,F12.4,A)') ' q (specific heat) = ', Q/mass, ' kJ/kg' write(*,'(A,F12.4,A)') ' du = ', dU/mass, ' kJ/kg' write(*,'(A,F12.4,A)') ' dh = ', dH/mass, ' kJ/kg' write(*,'(A,ES14.6,A)') ' ds = ', dS/mass, ' kJ/(kg.K)' end if write(*,*) ! ── Sensitivity sweep: n from 0.2 to 3.0 ─────────────────── n_sweep = 40 write(*,'(A)') '--- SENSITIVITY: POLYTROPIC INDEX SWEEP ---------------------' write(*,'(A)') ' n T2[K] W[kJ] Q[kJ]' write(*,'(A)') ' -----------------------------------------------------------' do i = 1, n_sweep n_sw = 0.2d0 + dble(i-1) * (3.0d0 - 0.2d0) / dble(n_sweep - 1) if (abs(n_sw - 1.0d0) < 0.02d0) then T2_sw = T1 V2_sw = V1 * (P1 / P2) W_sw = P1 * V1 * log(V2_sw / V1) else V2_sw = V1 * (P1/P2) ** (1.0d0/n_sw) T2_sw = T1 * (P2/P1) ** ((n_sw-1.0d0)/n_sw) W_sw = (P2*V2_sw - P1*V1) / (1.0d0 - n_sw) end if Q_sw = mass * cv * (T2_sw - T1) / 1000.0d0 + W_sw write(*,'(F8.4,4X,F12.2,4X,F12.4,4X,F12.4)') n_sw, T2_sw, W_sw, Q_sw end do write(*,*) write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------' write(*,'(A)') ' Polytropic relation: P*V^n = constant.' write(*,'(A)') ' V2 = V1*(P1/P2)^(1/n).' write(*,'(A)') ' T2 = T1*(P2/P1)^((n-1)/n).' write(*,'(A)') ' W = (P2V2 - P1V1)/(1-n) for n != 1.' write(*,'(A)') ' W = P1*V1*ln(V2/V1) for n = 1 (isothermal).' write(*,'(A)') ' Q = dU + W (first law, closed system).' write(*,'(A)') ' dS = m*[cv*ln(T2/T1) + R*ln(V2/V1)] (ideal gas).' write(*,'(A)') ' Special cases: n=0 isobaric, n=1 isothermal,' write(*,'(A)') ' n=gamma isentropic, n->inf isochoric.' end program isentropic_polytropic