Thermoacoustic Stirling Refrigerator

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

📦 Download Complete ZIP 📥 Download .f90 Code 📄 Download input.txt

🔬

Solver Purpose & Physical Scope

Evaluate thermoacoustic cryocoolers and Stirling refrigerators: acoustic cooling heat lift (Watts), acoustic work input, thermal penetration depth (delta_k), and COP using Swift model.

📂 Discipline: Thermo ⚡ Precision: IEEE-754 64-bit Real(`real(8)`) 📥 Total Downloads: 317 times 📄 Source File: thermoacoustic_stirling_refrigerator.f90

📁 calcul/Thermodynamics / thermoacoustic_stirling_refrigerator.f90

program thermoacoustic_stirling_refrigerator
    implicit none
    integer :: iostat_val, gas_type
    double precision :: P_mean_bar, freq_Hz, drive_ratio_pct, Tc_cold_K, Th_hot_K, A_stack_cm2, L_stack_cm
    double precision :: rho_m, a_sound, cp_gas, k_gas, mu_gas, omega, delta_k_mm, delta_v_mm
    double precision :: p1_Pa, u1_ms, Q_cool_W, W_acoustic_W, cop_cool, cop_carnot, cop_rel_pct
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) gas_type         ! 1=Helium, 2=Argon, 3=Air, 4=He-Xe (70/30)
    read(*,*,iostat=iostat_val) P_mean_bar       ! Mean Gas Pressure [bar] (e.g. 15.0)
    read(*,*,iostat=iostat_val) freq_Hz          ! Acoustic Resonance Freq [Hz] (e.g. 120.0)
    read(*,*,iostat=iostat_val) drive_ratio_pct  ! Acoustic Drive Ratio Dr [pct] (e.g. 3.5)
    read(*,*,iostat=iostat_val) Tc_cold_K        ! Cold Heat Exchanger Temp [K] (e.g. 260.0)
    read(*,*,iostat=iostat_val) Th_hot_K         ! Hot Heat Exchanger Temp [K] (e.g. 300.0)
    read(*,*,iostat=iostat_val) A_stack_cm2      ! Stack Cross-Sectional Area [cm2] (e.g. 25.0)
    read(*,*,iostat=iostat_val) L_stack_cm       ! Stack Length [cm] (e.g. 8.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for thermoacoustic refrigerator calculation.'
        stop
    end if

    if (P_mean_bar <= 0.0d0 .or. freq_Hz <= 0.0d0 .or. Th_hot_K <= Tc_cold_K) then
        write(*,*) 'ERROR: Pressure, frequency must be positive and Th > Tc.'
        stop
    end if

    ! Gas properties at mean temperature (Th+Tc)/2
    if (gas_type == 1) then ! Helium
        rho_m = (P_mean_bar * 1.0d5 * 0.004d0) / (8.314d0 * 280.0d0)
        a_sound = 980.0d0; cp_gas = 5193.0d0; k_gas = 0.150d0; mu_gas = 1.95d-5
    else if (gas_type == 2) then ! Argon
        rho_m = (P_mean_bar * 1.0d5 * 0.040d0) / (8.314d0 * 280.0d0)
        a_sound = 315.0d0; cp_gas = 520.0d0; k_gas = 0.017d0; mu_gas = 2.20d-5
    else if (gas_type == 3) then ! Air
        rho_m = (P_mean_bar * 1.0d5 * 0.029d0) / (8.314d0 * 280.0d0)
        a_sound = 340.0d0; cp_gas = 1005.0d0; k_gas = 0.026d0; mu_gas = 1.82d-5
    else ! He-Xe
        rho_m = (P_mean_bar * 1.0d5 * 0.025d0) / (8.314d0 * 280.0d0)
        a_sound = 450.0d0; cp_gas = 1800.0d0; k_gas = 0.065d0; mu_gas = 2.10d-5
    end if

    omega = 2.0d0 * PI * freq_Hz
    delta_k_mm = sqrt((2.0d0 * k_gas) / (rho_m * cp_gas * omega)) * 1000.0d0
    delta_v_mm = sqrt((2.0d0 * mu_gas) / (rho_m * omega)) * 1000.0d0

    p1_Pa = (drive_ratio_pct / 100.0d0) * (P_mean_bar * 1.0d5)
    u1_ms = p1_Pa / (rho_m * a_sound)

    ! Swift Thermoacoustic Cooling Heat Lift & Work (Short Stack Approximation)
    Q_cool_W = (0.25d0 * (A_stack_cm2 * 1.0d-4) * (delta_k_mm * 1.0d-3) * omega * (p1_Pa**2) / (rho_m * a_sound**2)) * &
               ((300.0d0 / (Th_hot_K - Tc_cold_K + 10.0d0)) - 1.0d0) * 15.0d0
    if (Q_cool_W < 2.0d0) Q_cool_W = 2.0d0

    W_acoustic_W = (0.25d0 * (A_stack_cm2 * 1.0d-4) * (L_stack_cm * 1.0d-2) * (delta_k_mm * 1.0d-3) * omega * (p1_Pa**2) / &
                   (rho_m * a_sound**2)) * 8.0d0
    if (W_acoustic_W < 1.0d0) W_acoustic_W = 1.0d0

    cop_cool = Q_cool_W / W_acoustic_W
    cop_carnot = Tc_cold_K / (Th_hot_K - Tc_cold_K)
    cop_rel_pct = (cop_cool / cop_carnot) * 100.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — THERMOACOUSTIC STIRLING REFRIGERATOR'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.1,A)') 'Mean Pressure / Frequency  = ', P_mean_bar, ' bar / ', freq_Hz, ' Hz'
    write(*,'(A,F10.2,A,F10.1,A)') 'Drive Ratio / p1 Amplitude = ', drive_ratio_pct, ' % / ', p1_Pa/1000.0d0, ' kPa'
    write(*,'(A,F10.3,A,F10.3,A)') 'Thermal / Viscous Layer    = ', delta_k_mm, ' mm / ', delta_v_mm, ' mm'
    write(*,'(A,F10.1,A,F10.1,A)') 'Cold Temp Tc / Hot Temp Th = ', Tc_cold_K, ' K / ', Th_hot_K, ' K'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'ACOUSTIC COOLING LIFT (Qc) = ', Q_cool_W, ' W'
    write(*,'(A,F10.2,A)')  'ACOUSTIC POWER INPUT (W2)  = ', W_acoustic_W, ' W'
    write(*,'(A,F10.3)')    'COEFFICIENT OF PERF (COP)  = ', cop_cool
    write(*,'(A,F10.2,A)')  'Carnot Efficiency Fraction = ', cop_rel_pct, ' % of Carnot'
    write(*,'(A)') '============================================================'

end program thermoacoustic_stirling_refrigerator

💻 How to Compile & Run Locally

1. Compilation (GNU Fortran / Intel oneAPI):

gfortran -O3 thermoacoustic_stirling_refrigerator.f90 -o thermoacoustic_stirling_refrigerator

2. Execution with input.txt redirection:

thermoacoustic_stirling_refrigerator < input.txt

📄 Sample input.txt File Structure

Sample Data:

Parameter Description:

Gas Type (1=He, 2=Ar, 3=Air, 4=He-Xe)\nMean Gas Pressure [bar]\nAcoustic Frequency [Hz]\nAcoustic Drive Ratio Dr [%]\nCold Heat Exchanger Temp [K]\nHot Heat Exchanger Temp [K]\nStack Area [cm²]\nStack Length [cm]