⚡ Fortran 90 / 2008 Double Precision
📥 280 Downloads

Carnot Cycle & Heat Pump Limit

Standalone, self-contained numerical routine. Verify algorithms, inspect boundary condition equations, or compile locally for batch parametric runs.

🔬

Solver Purpose & Physical Scope

Solves the theoretical limits of thermodynamic cycles operating between hot and cold temperature reservoirs. Calculates Carnot efficiency for heat engines, Carnot COP for refrigerators and heat pumps, and estimates real cycle performance metrics (net work, heat exchange, actual efficiency/COP) using a user-defined fraction of Carnot (real cycle factor).

📂 Discipline: Thermo Precision: IEEE-754 64-bit Real(`real(8)`) 📥 Total Downloads: 280 times 📄 Source File: carnot_heatpump.f90
📁 calcul/Thermodynamics / carnot_heatpump.f90
program carnot_heatpump
    implicit none
    integer :: cycle_type, iostat_val, i, n_sweep
    double precision :: T_hot, T_cold, Q_desired, real_frac
    double precision :: eta_carnot, COP_ref, COP_hp
    double precision :: W_cycle, Q_hot, Q_cold
    double precision :: COP_ref_real, COP_hp_real, eta_real
    double precision :: W_real, Q_hot_real, Q_cold_real
    double precision :: T_sw, eta_sw, COP_ref_sw, COP_hp_sw
    character(len=40) :: cycle_name

    ! ── Read inputs ─────────────────────────────────────────────
    read(*,*,iostat=iostat_val) T_hot
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid T_hot input.'
        stop
    end if
    read(*,*,iostat=iostat_val) T_cold
    read(*,*,iostat=iostat_val) Q_desired
    read(*,*,iostat=iostat_val) cycle_type
    read(*,*,iostat=iostat_val) real_frac
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all inputs.'
        stop
    end if

    ! ── Input validation ────────────────────────────────────────
    if (T_hot <= 0.0d0) then
        write(*,*) 'ERROR: T_hot must be positive (K).'
        stop
    end if
    if (T_cold <= 0.0d0) then
        write(*,*) 'ERROR: T_cold must be positive (K).'
        stop
    end if
    if (T_hot <= T_cold) then
        write(*,*) 'ERROR: T_hot must exceed T_cold.'
        stop
    end if
    if (Q_desired <= 0.0d0) Q_desired = 100.0d0
    if (cycle_type < 1 .or. cycle_type > 3) cycle_type = 1
    if (real_frac <= 0.0d0 .or. real_frac > 1.0d0) real_frac = 0.50d0

    if (cycle_type == 1) then
        cycle_name = 'Carnot Heat Engine'
    else if (cycle_type == 2) then
        cycle_name = 'Refrigerator (Reversed Carnot)'
    else
        cycle_name = 'Heat Pump (Reversed Carnot)'
    end if

    ! ── Carnot limits ───────────────────────────────────────────
    eta_carnot = 1.0d0 - T_cold / T_hot
    COP_ref = T_cold / (T_hot - T_cold)
    COP_hp  = T_hot  / (T_hot - T_cold)

    ! ── Energy balance based on cycle type ──────────────────────
    if (cycle_type == 1) then
        ! Heat engine: Q_desired = Q_hot (heat input)
        Q_hot = Q_desired
        W_cycle = Q_hot * eta_carnot
        Q_cold = Q_hot - W_cycle
    else if (cycle_type == 2) then
        ! Refrigerator: Q_desired = Q_cold (cooling load)
        Q_cold = Q_desired
        W_cycle = Q_cold / max(COP_ref, 1.0d-10)
        Q_hot = Q_cold + W_cycle
    else
        ! Heat pump: Q_desired = Q_hot (heating load)
        Q_hot = Q_desired
        W_cycle = Q_hot / max(COP_hp, 1.0d-10)
        Q_cold = Q_hot - W_cycle
    end if

    ! ── Real cycle estimates ────────────────────────────────────
    eta_real = eta_carnot * real_frac
    COP_ref_real = COP_ref * real_frac
    COP_hp_real  = COP_hp  * real_frac

    if (cycle_type == 1) then
        W_real = Q_desired * eta_real
        Q_hot_real = Q_desired
        Q_cold_real = Q_hot_real - W_real
    else if (cycle_type == 2) then
        W_real = Q_desired / max(COP_ref_real, 1.0d-10)
        Q_cold_real = Q_desired
        Q_hot_real = Q_cold_real + W_real
    else
        W_real = Q_desired / max(COP_hp_real, 1.0d-10)
        Q_hot_real = Q_desired
        Q_cold_real = Q_hot_real - W_real
    end if

    ! ── Output ──────────────────────────────────────────────────
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   CARNOT, REVERSED CARNOT & HEAT PUMP'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,F12.2,A)')  '  Hot Reservoir T_hot       = ', T_hot, ' K'
    write(*,'(A,F12.2,A)')  '  Hot Reservoir T_hot       = ', T_hot-273.15d0, ' C'
    write(*,'(A,F12.2,A)')  '  Cold Reservoir T_cold     = ', T_cold, ' K'
    write(*,'(A,F12.2,A)')  '  Cold Reservoir T_cold     = ', T_cold-273.15d0, ' C'
    write(*,'(A,F12.2,A)')  '  Q_desired                 = ', Q_desired, ' kW'
    write(*,'(A,A)')        '  Cycle Type                = ', trim(cycle_name)
    write(*,'(A,F10.4)')    '  Real Cycle Factor         = ', real_frac
    write(*,*)
    write(*,'(A)') '--- CARNOT LIMITS -------------------------------------------'
    write(*,'(A,F10.6)')    '  Carnot Efficiency         = ', eta_carnot
    write(*,'(A,F10.2,A)')  '  Carnot Efficiency         = ', eta_carnot*100.0d0, ' percent'
    write(*,'(A,F10.4)')    '  COP Refrigeration (max)   = ', COP_ref
    write(*,'(A,F10.4)')    '  COP Heat Pump (max)       = ', COP_hp
    write(*,'(A)')          '  COP_hp = COP_ref + 1      (verified)'
    write(*,*)
    write(*,'(A)') '--- ENERGY BALANCE (IDEAL) ----------------------------------'
    write(*,'(A,F12.4,A)')  '  Q_hot                     = ', Q_hot, ' kW'
    write(*,'(A,F12.4,A)')  '  Q_cold                    = ', Q_cold, ' kW'
    write(*,'(A,F12.4,A)')  '  W_cycle                   = ', W_cycle, ' kW'
    write(*,*)
    write(*,'(A)') '--- REAL CYCLE ESTIMATES ------------------------------------'
    write(*,'(A,F10.4)')    '  Real Efficiency           = ', eta_real
    write(*,'(A,F10.2,A)')  '  Real Efficiency           = ', eta_real*100.0d0, ' percent'
    write(*,'(A,F10.4)')    '  Real COP Refrigeration    = ', COP_ref_real
    write(*,'(A,F10.4)')    '  Real COP Heat Pump        = ', COP_hp_real
    write(*,'(A,F12.4,A)')  '  Real W_cycle              = ', W_real, ' kW'
    write(*,'(A,F12.4,A)')  '  Real Q_hot                = ', Q_hot_real, ' kW'
    write(*,'(A,F12.4,A)')  '  Real Q_cold               = ', Q_cold_real, ' kW'
    write(*,*)

    ! ── Sensitivity sweep: T_cold from 200 K to T_hot-5 ────────
    n_sweep = 40
    write(*,'(A)') '--- SENSITIVITY: PERFORMANCE VS T_COLD ----------------------'
    write(*,'(A)') '  T_cold[K]     eta_Carnot    COP_ref       COP_hp'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, n_sweep
        T_sw = 200.0d0 + dble(i-1) * (T_hot - 5.0d0 - 200.0d0) / dble(n_sweep - 1)
        if (T_sw >= T_hot) T_sw = T_hot - 1.0d0
        eta_sw = 1.0d0 - T_sw / T_hot
        COP_ref_sw = T_sw / max(T_hot - T_sw, 1.0d-10)
        COP_hp_sw  = T_hot / max(T_hot - T_sw, 1.0d-10)
        write(*,'(F10.2,4X,F10.6,4X,F10.4,4X,F10.4)') T_sw, eta_sw, COP_ref_sw, COP_hp_sw
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Carnot efficiency: eta = 1 - T_cold/T_hot'
    write(*,'(A)') '  COP_refrigeration = T_cold / (T_hot - T_cold)'
    write(*,'(A)') '  COP_heat_pump = T_hot / (T_hot - T_cold)'
    write(*,'(A)') '  COP_hp = COP_ref + 1'
    write(*,'(A)') '  Real cycle: multiply ideal COP by factor (0.4-0.6 typical)'
    write(*,'(A)') '  Energy balance: Q_hot = Q_cold + W'

end program carnot_heatpump


💻 How to Compile & Run Locally

1. Compilation (GNU Fortran / Intel oneAPI):

gfortran -O3 carnot_heatpump.f90 -o carnot_heatpump

2. Execution with input.txt redirection:

carnot_heatpump < input.txt

📄 Sample input.txt File Structure

Sample Data:
773.15
303.15
500000.0
1
0.50
Parameter Description:
Hot reservoir temperature TH [K]
Cold reservoir temperature TL [K]
Desired heat or power Q [kW]
Cycle type (1=Power Cycle, 2=Refrigeration, 3=Heat Pump)
Real cycle factor (fraction of Carnot)