🔬
Solver Purpose & Physical Scope
Model ideal Stirling and Ericsson gas cycles operating with perfect regeneration. Computes cycle state points (pressure, volume, temperature), net boundary work, thermal efficiency with and without regeneration, maximum and minimum cycle pressures, and equivalent refrigeration/heating COP limits. Generates P-V and T-s diagram coordinate curves.
📂 Discipline: Thermo
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 404 times
📄 Source File:
stirling_ericsson.f90
📁 calcul/Thermodynamics /
stirling_ericsson.f90
program stirling_ericsson
implicit none
integer :: cycle_type, i, iostat_val, n_pts
double precision :: TH, TL, Pmax, Pmin, Vmax, Vmin, R_gas, cv, cp, gamma
double precision :: Mmol, n_mol
double precision :: Wnet, Qin, Qout, Qregen, eta_th, eta_carnot
double precision :: r_comp, r_press, COP_heat, COP_ref
double precision :: W12, W23, W34, W41, Q12, Q23, Q34, Q41
double precision :: P1, P2, P3, P4, V1, V2, V3, V4, T1, T2, T3, T4
double precision :: theta, P_i, V_i, T_i
double precision, parameter :: PI = 3.141592653589793d0
double precision, parameter :: Ru = 8.314462d0
character(len=40) :: cycle_name
read(*,*,iostat=iostat_val) cycle_type
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid cycle type input.'
stop
end if
read(*,*,iostat=iostat_val) TH
read(*,*,iostat=iostat_val) TL
read(*,*,iostat=iostat_val) Vmin
read(*,*,iostat=iostat_val) Vmax
read(*,*,iostat=iostat_val) Mmol
read(*,*,iostat=iostat_val) gamma
read(*,*,iostat=iostat_val) n_mol
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all cycle inputs.'
stop
end if
if (TH <= TL .or. TH <= 0.0d0 .or. TL <= 0.0d0) then
write(*,*) 'ERROR: TH must be greater than TL, both positive.'
stop
end if
if (Vmin <= 0.0d0 .or. Vmax <= Vmin) then
write(*,*) 'ERROR: Volumes must be positive, Vmax > Vmin.'
stop
end if
if (Mmol <= 0.0d0) Mmol = 4.003d0 ! helium default
if (gamma <= 1.0d0) gamma = 1.667d0
if (n_mol <= 0.0d0) n_mol = 1.0d0
R_gas = Ru ! per mole
cv = R_gas / (gamma - 1.0d0)
cp = gamma * cv
r_comp = Vmax / Vmin
eta_carnot = 1.0d0 - TL / TH
select case(cycle_type)
case(1)
! ====== STIRLING CYCLE ======
! 1-2: Isothermal compression at TL (V1=Vmax → V2=Vmin)
! 2-3: Isochoric heating at Vmin (TL → TH)
! 3-4: Isothermal expansion at TH (V3=Vmin → V4=Vmax)
! 4-1: Isochoric cooling at Vmax (TH → TL)
cycle_name = 'Stirling Cycle'
T1 = TL; V1 = Vmax; P1 = n_mol*R_gas*T1/V1
T2 = TL; V2 = Vmin; P2 = n_mol*R_gas*T2/V2
T3 = TH; V3 = Vmin; P3 = n_mol*R_gas*T3/V3
T4 = TH; V4 = Vmax; P4 = n_mol*R_gas*T4/V4
! Work (per mole * n_mol)
W12 = n_mol * R_gas * TL * log(V2/V1) ! negative (compression)
W23 = 0.0d0 ! isochoric
W34 = n_mol * R_gas * TH * log(V4/V3) ! positive (expansion)
W41 = 0.0d0 ! isochoric
! Heat
Q12 = W12 ! isothermal: Q = W
Q23 = n_mol * cv * (TH - TL) ! isochoric heating
Q34 = W34 ! isothermal: Q = W
Q41 = n_mol * cv * (TL - TH) ! isochoric cooling
! With perfect regenerator: Q23 recovered from Q41
Qregen = n_mol * cv * (TH - TL)
case(2)
! ====== ERICSSON CYCLE ======
! 1-2: Isothermal compression at TL (P1=Pmin → P2=Pmax)
! 2-3: Isobaric heating at Pmax (TL → TH)
! 3-4: Isothermal expansion at TH (P3=Pmax → P4=Pmin)
! 4-1: Isobaric cooling at Pmin (TH → TL)
cycle_name = 'Ericsson Cycle'
Pmin = n_mol * R_gas * TL / Vmax
Pmax = n_mol * R_gas * TL / Vmin
r_press = Pmax / Pmin
T1 = TL; P1 = Pmin; V1 = n_mol*R_gas*T1/P1
T2 = TL; P2 = Pmax; V2 = n_mol*R_gas*T2/P2
T3 = TH; P3 = Pmax; V3 = n_mol*R_gas*T3/P3
T4 = TH; P4 = Pmin; V4 = n_mol*R_gas*T4/P4
! Work
W12 = n_mol * R_gas * TL * log(V2/V1) ! negative
W23 = P3 * (V3 - V2) ! isobaric expansion work (by gas during heating)
W34 = n_mol * R_gas * TH * log(V4/V3) ! positive
W41 = P1 * (V1 - V4) ! isobaric compression work (negative)
! Heat
Q12 = W12 ! isothermal
Q23 = n_mol * cp * (TH - TL) ! isobaric heating
Q34 = W34 ! isothermal
Q41 = n_mol * cp * (TL - TH) ! isobaric cooling
Qregen = n_mol * cp * (TH - TL)
case default
write(*,*) 'ERROR: Cycle type must be 1 (Stirling) or 2 (Ericsson).'
stop
end select
Wnet = W12 + W23 + W34 + W41
Qin = Q34 ! heat added during hot isothermal
if (cycle_type == 2) then
! Without regenerator, Qin includes isobaric heating too
! With perfect regenerator, only isothermal heat input counts
end if
Qout = abs(Q12)
! With perfect regenerator: eta = eta_Carnot
! Without regenerator:
if (cycle_type == 1) then
! Without regen: Qin_total = Q34 + Q23
eta_th = Wnet / (Q34 + Q23)
else
eta_th = Wnet / (Q34 + Q23)
end if
! COP as heat pump
COP_heat = (Q34 + Q23) / max(abs(Wnet), 1.0d-30)
! COP as refrigerator
COP_ref = Qout / max(abs(Wnet), 1.0d-30)
Pmax = max(P1, max(P2, max(P3, P4)))
Pmin = min(P1, min(P2, min(P3, P4)))
write(*,'(A)') '============================================================'
write(*,'(A)') ' STIRLING & ERICSSON CYCLES ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUTS --------------------------------------------------'
write(*,'(A,A)') ' Cycle Type = ', trim(cycle_name)
write(*,'(A,ES12.4,A)') ' Hot Temperature TH = ', TH, ' K'
write(*,'(A,ES12.4,A)') ' Cold Temperature TL = ', TL, ' K'
write(*,'(A,ES12.4,A)') ' Min Volume = ', Vmin, ' m3'
write(*,'(A,ES12.4,A)') ' Max Volume = ', Vmax, ' m3'
write(*,'(A,ES12.4)') ' Compression Ratio = ', r_comp
write(*,'(A,ES12.4)') ' Gamma = ', gamma
write(*,'(A,ES12.4,A)') ' Molar Mass = ', Mmol, ' g/mol'
write(*,'(A,ES12.4)') ' Moles of Gas = ', n_mol
write(*,*)
write(*,'(A)') '--- STATE POINTS --------------------------------------------'
write(*,'(A)') ' Point T[K] P[Pa] V[m3]'
write(*,'(A)') ' -------------------------------------------'
write(*,'(A,F12.2,2X,ES12.4,2X,ES12.4)') ' 1 ', T1, P1, V1
write(*,'(A,F12.2,2X,ES12.4,2X,ES12.4)') ' 2 ', T2, P2, V2
write(*,'(A,F12.2,2X,ES12.4,2X,ES12.4)') ' 3 ', T3, P3, V3
write(*,'(A,F12.2,2X,ES12.4,2X,ES12.4)') ' 4 ', T4, P4, V4
write(*,*)
write(*,'(A)') '--- PROCESS ENERGY ------------------------------------------'
write(*,'(A,ES12.4,A)') ' W12 (1→2) = ', W12, ' J'
write(*,'(A,ES12.4,A)') ' Q12 (1→2) = ', Q12, ' J'
write(*,'(A,ES12.4,A)') ' W23 (2→3) = ', W23, ' J'
write(*,'(A,ES12.4,A)') ' Q23 (2→3) = ', Q23, ' J'
write(*,'(A,ES12.4,A)') ' W34 (3→4) = ', W34, ' J'
write(*,'(A,ES12.4,A)') ' Q34 (3→4) = ', Q34, ' J'
write(*,'(A,ES12.4,A)') ' W41 (4→1) = ', W41, ' J'
write(*,'(A,ES12.4,A)') ' Q41 (4→1) = ', Q41, ' J'
write(*,*)
write(*,'(A)') '--- CYCLE PERFORMANCE ---------------------------------------'
write(*,'(A,ES12.4,A)') ' Net Work Wnet = ', Wnet, ' J'
write(*,'(A,ES12.4,A)') ' Heat Input (no regen) = ', Q34+Q23, ' J'
write(*,'(A,ES12.4,A)') ' Heat Input (with regen) = ', Q34, ' J'
write(*,'(A,ES12.4,A)') ' Heat Rejected = ', abs(Q12)+abs(Q41), ' J'
write(*,'(A,ES12.4,A)') ' Regenerator Heat = ', Qregen, ' J'
write(*,'(A,F12.4)') ' Thermal Efficiency (no regen) = ', eta_th
write(*,'(A,F12.4)') ' Thermal Efficiency (with regen)= ', eta_carnot
write(*,'(A,F12.4)') ' Carnot Efficiency = ', eta_carnot
write(*,'(A,F12.4)') ' COP (heat pump) = ', COP_heat
write(*,'(A,F12.4)') ' COP (refrigerator) = ', COP_ref
write(*,'(A,ES12.4,A)') ' Max Pressure = ', Pmax, ' Pa'
write(*,'(A,ES12.4,A)') ' Min Pressure = ', Pmin, ' Pa'
write(*,*)
! P-V diagram data
n_pts = 50
write(*,'(A)') '--- PV DIAGRAM DATA -----------------------------------------'
write(*,'(A)') ' V[m3] P[Pa] process'
write(*,'(A)') ' -----------------------------------------------------------'
if (cycle_type == 1) then
! 1→2: isothermal TL, V1→V2
do i = 0, n_pts
V_i = V1 + (V2-V1)*dble(i)/dble(n_pts)
P_i = n_mol*R_gas*TL/V_i
write(*,'(ES12.4,2X,ES12.4,2X,A)') V_i, P_i, '1-2 iso-T(TL)'
end do
! 2→3: isochoric, V=Vmin
do i = 0, n_pts
T_i = TL + (TH-TL)*dble(i)/dble(n_pts)
P_i = n_mol*R_gas*T_i/V2
write(*,'(ES12.4,2X,ES12.4,2X,A)') V2, P_i, '2-3 iso-V'
end do
! 3→4: isothermal TH, V3→V4
do i = 0, n_pts
V_i = V3 + (V4-V3)*dble(i)/dble(n_pts)
P_i = n_mol*R_gas*TH/V_i
write(*,'(ES12.4,2X,ES12.4,2X,A)') V_i, P_i, '3-4 iso-T(TH)'
end do
! 4→1: isochoric, V=Vmax
do i = 0, n_pts
T_i = TH + (TL-TH)*dble(i)/dble(n_pts)
P_i = n_mol*R_gas*T_i/V4
write(*,'(ES12.4,2X,ES12.4,2X,A)') V4, P_i, '4-1 iso-V'
end do
else
! Ericsson
! 1→2: isothermal TL
do i = 0, n_pts
V_i = V1 + (V2-V1)*dble(i)/dble(n_pts)
P_i = n_mol*R_gas*TL/V_i
write(*,'(ES12.4,2X,ES12.4,2X,A)') V_i, P_i, '1-2 iso-T(TL)'
end do
! 2→3: isobaric Pmax
do i = 0, n_pts
T_i = TL + (TH-TL)*dble(i)/dble(n_pts)
V_i = n_mol*R_gas*T_i/P3
write(*,'(ES12.4,2X,ES12.4,2X,A)') V_i, P3, '2-3 iso-P'
end do
! 3→4: isothermal TH
do i = 0, n_pts
V_i = V3 + (V4-V3)*dble(i)/dble(n_pts)
P_i = n_mol*R_gas*TH/V_i
write(*,'(ES12.4,2X,ES12.4,2X,A)') V_i, P_i, '3-4 iso-T(TH)'
end do
! 4→1: isobaric Pmin
do i = 0, n_pts
T_i = TH + (TL-TH)*dble(i)/dble(n_pts)
V_i = n_mol*R_gas*T_i/P1
write(*,'(ES12.4,2X,ES12.4,2X,A)') V_i, P1, '4-1 iso-P'
end do
end if
write(*,*)
! T-s diagram data
write(*,'(A)') '--- TS DIAGRAM DATA -----------------------------------------'
write(*,'(A)') ' s_rel[J/K] T[K] process'
write(*,'(A)') ' -----------------------------------------------------------'
call write_ts_data(cycle_type, T1, T2, T3, T4, V1, V2, V3, V4, &
P1, P2, P3, P4, n_mol, R_gas, cv, cp, n_pts)
write(*,*)
! Efficiency vs temperature ratio sweep
write(*,'(A)') '--- EFFICIENCY VS TL/TH SWEEP -------------------------------'
write(*,'(A)') ' TL/TH eta_Carnot eta_no_regen'
write(*,'(A)') ' -------------------------------------------'
do i = 1, 40
theta = 0.1d0 + 0.85d0*dble(i-1)/39.0d0
write(*,'(F10.4,2X,F10.5,2X,F10.5)') theta, 1.0d0-theta, &
(1.0d0-theta)*log(r_comp) / &
(log(r_comp) + cv/R_gas*(1.0d0/theta - 1.0d0)/(1.0d0))
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' Stirling: iso-T compression/expansion + iso-V regeneration.'
write(*,'(A)') ' Ericsson: iso-T compression/expansion + iso-P regeneration.'
write(*,'(A)') ' With perfect regen: eta = eta_Carnot = 1 - TL/TH.'
contains
subroutine write_ts_data(ctype, t1i, t2i, t3i, t4i, v1i, v2i, v3i, v4i, &
p1i, p2i, p3i, p4i, nm, Rg, cvi, cpi, npts)
implicit none
integer, intent(in) :: ctype, npts
double precision, intent(in) :: t1i,t2i,t3i,t4i,v1i,v2i,v3i,v4i
double precision, intent(in) :: p1i,p2i,p3i,p4i,nm,Rg,cvi,cpi
double precision :: s_acc, ds_i, T_loc, frac
integer :: j
s_acc = 0.0d0
if (ctype == 1) then
! 1→2: isothermal at TL
do j = 0, npts
frac = dble(j)/dble(npts)
ds_i = nm*Rg*log(v1i + (v2i-v1i)*frac) - nm*Rg*log(v1i)
write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, t1i, '1-2'
end do
s_acc = s_acc + nm*Rg*log(v2i/v1i)
! 2→3: isochoric
do j = 0, npts
frac = dble(j)/dble(npts)
T_loc = t2i + (t3i-t2i)*frac
ds_i = nm*cvi*log(T_loc/t2i)
write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, T_loc, '2-3'
end do
s_acc = s_acc + nm*cvi*log(t3i/t2i)
! 3→4: isothermal at TH
do j = 0, npts
frac = dble(j)/dble(npts)
ds_i = nm*Rg*log(v3i + (v4i-v3i)*frac) - nm*Rg*log(v3i)
write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, t3i, '3-4'
end do
s_acc = s_acc + nm*Rg*log(v4i/v3i)
! 4→1: isochoric
do j = 0, npts
frac = dble(j)/dble(npts)
T_loc = t4i + (t1i-t4i)*frac
ds_i = nm*cvi*log(T_loc/t4i)
write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, T_loc, '4-1'
end do
else
! Ericsson
! 1→2: isothermal
do j = 0, npts
frac = dble(j)/dble(npts)
ds_i = nm*Rg*log(v1i + (v2i-v1i)*frac) - nm*Rg*log(v1i)
write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, t1i, '1-2'
end do
s_acc = s_acc + nm*Rg*log(v2i/v1i)
! 2→3: isobaric
do j = 0, npts
frac = dble(j)/dble(npts)
T_loc = t2i + (t3i-t2i)*frac
ds_i = nm*cpi*log(T_loc/t2i)
write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, T_loc, '2-3'
end do
s_acc = s_acc + nm*cpi*log(t3i/t2i)
! 3→4: isothermal
do j = 0, npts
frac = dble(j)/dble(npts)
ds_i = nm*Rg*log(v3i + (v4i-v3i)*frac) - nm*Rg*log(v3i)
write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, t3i, '3-4'
end do
s_acc = s_acc + nm*Rg*log(v4i/v3i)
! 4→1: isobaric
do j = 0, npts
frac = dble(j)/dble(npts)
T_loc = t4i + (t1i-t4i)*frac
ds_i = nm*cpi*log(T_loc/t4i)
write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, T_loc, '4-1'
end do
end if
end subroutine write_ts_data
end program stirling_ericsson
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 stirling_ericsson.f90 -o stirling_ericsson
2. Execution with input.txt redirection:
stirling_ericsson < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 900.0 300.0 0.0005 0.002 4.003 1.667 1.0
Parameter Description:
Cycle type (1=Stirling, 2=Ericsson) Hot source temperature TH [K] Cold sink temperature TL [K] Minimum gas volume Vmin [m3] Maximum gas volume Vmax [m3] Molar mass of gas [g/mol] Specific heat ratio gamma (cp/cv) Amount of gas [moles]