🔬
Solver Purpose & Physical Scope
Analyzes ideal gas expansion and compression processes under polytropic and isentropic constraints. Solves final thermodynamic states (pressure, temperature, volume) and computes boundary work, heat transfer, internal energy change ($\Delta U$), enthalpy change ($\Delta H$), and entropy change ($\Delta S$) for standard thermodynamic processes (isochoric, isobaric, isothermal, isentropic, polytropic).
📂 Discipline: Thermo
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 509 times
📄 Source File:
isentropic_polytropic.f90
📁 calcul/Thermodynamics /
isentropic_polytropic.f90
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 isentropic_polytropic.f90 -o isentropic_polytropic
2. Execution with input.txt redirection:
isentropic_polytropic < input.txt
📄 Sample input.txt File Structure
Sample Data:
100.0 800.0 300.0 0.5 1.4 287.0 1.4 1.0
Parameter Description:
Initial pressure P1 [kPa] Final pressure P2 [kPa] Initial temperature T1 [K] Initial volume V1 [m3] Polytropic index n Specific gas constant R [J/kg-K] Specific heat ratio gamma (cp/cv) Gas mass [kg]