Two-Phase Closed Thermosiphon & Heat Pipe Sizer

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

Size two-phase closed thermosiphons (TPCT) and heat pipes: flooding/entrainment limit, boiling burnout limit, effective thermal conductivity (W/m·K), and temperature drop.

📂 Discipline: Heat-exchangers ⚡ Precision: IEEE-754 64-bit Real(`real(8)`) 📥 Total Downloads: 400 times 📄 Source File: two_phase_closed_thermosiphon_heat_pipe.f90

📁 calcul/HeatExchangers / two_phase_closed_thermosiphon_heat_pipe.f90

program two_phase_closed_thermosiphon_heat_pipe
    implicit none
    integer :: iostat_val, fluid_type, pipe_material
    double precision :: Q_heat_W, Tv_sat_C, D_outer_mm, D_inner_mm, Le_mm, La_mm, Lc_mm, fill_ratio_pct
    double precision :: rho_L, rho_v, cp_L, mu_L, mu_v, k_L, hfg_kJkg, sigma_Nm, Pv_bar
    double precision :: D_i_m, D_o_m, L_e_m, L_a_m, L_c_m, L_eff_m, A_v_m2, A_e_m2, A_c_m2
    double precision :: Q_flood_W, Q_boiling_W, Q_sonic_W, Q_max_allow_W, delta_T_pipe_C, k_eff_WmK
    character(len=32) :: limit_type, status_str
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: G_ACC = 9.80665d0

    ! Read inputs
    read(*,*,iostat=iostat_val) fluid_type      ! 1=Water, 2=Ethanol, 3=Methanol, 4=Acetone, 5=Ammonia
    read(*,*,iostat=iostat_val) pipe_material   ! 1=Copper, 2=Aluminum, 3=Stainless Steel
    read(*,*,iostat=iostat_val) Q_heat_W        ! Transferred Thermal Load [W] (e.g. 500.0)
    read(*,*,iostat=iostat_val) Tv_sat_C        ! Operating Saturation Temp [deg C] (e.g. 70.0)
    read(*,*,iostat=iostat_val) D_outer_mm      ! Outer Pipe Diameter [mm] (e.g. 25.0)
    read(*,*,iostat=iostat_val) D_inner_mm      ! Inner Vapor Core Diameter [mm] (e.g. 22.0)
    read(*,*,iostat=iostat_val) Le_mm           ! Evaporator Length [mm] (e.g. 200.0)
    read(*,*,iostat=iostat_val) La_mm           ! Adiabatic Transport Length [mm] (e.g. 300.0)
    read(*,*,iostat=iostat_val) Lc_mm           ! Condenser Length [mm] (e.g. 200.0)
    read(*,*,iostat=iostat_val) fill_ratio_pct  ! Working Fluid Liquid Fill [% of Le] (e.g. 45.0)

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

    if (Q_heat_W <= 0.0d0 .or. D_inner_mm <= 0.0d0 .or. Le_mm <= 0.0d0 .or. Lc_mm <= 0.0d0) then
        write(*,*) 'ERROR: Thermal power, diameter, and section lengths must be positive.'
        stop
    end if

    ! Fluid thermophysical properties
    if (fluid_type == 1) then ! Water
        rho_L = 978.0d0; rho_v = 0.20d0; cp_L = 4190.0d0; mu_L = 0.00040d0
        k_L = 0.66d0; hfg_kJkg = 2330.0d0; sigma_Nm = 0.064d0; Pv_bar = 0.31d0
    else if (fluid_type == 2) then ! Ethanol
        rho_L = 750.0d0; rho_v = 1.45d0; cp_L = 2650.0d0; mu_L = 0.00055d0
        k_L = 0.16d0; hfg_kJkg = 880.0d0;  sigma_Nm = 0.018d0; Pv_bar = 0.72d0
    else if (fluid_type == 3) then ! Methanol
        rho_L = 740.0d0; rho_v = 1.30d0; cp_L = 2700.0d0; mu_L = 0.00035d0
        k_L = 0.19d0; hfg_kJkg = 1100.0d0; sigma_Nm = 0.019d0; Pv_bar = 1.25d0
    else if (fluid_type == 4) then ! Acetone
        rho_L = 740.0d0; rho_v = 2.10d0; cp_L = 2300.0d0; mu_L = 0.00025d0
        k_L = 0.16d0; hfg_kJkg = 520.0d0;  sigma_Nm = 0.019d0; Pv_bar = 1.60d0
    else ! Ammonia
        rho_L = 580.0d0; rho_v = 18.5d0; cp_L = 4900.0d0; mu_L = 0.00012d0
        k_L = 0.48d0; hfg_kJkg = 1180.0d0; sigma_Nm = 0.018d0; Pv_bar = 33.0d0
    end if

    D_i_m = D_inner_mm / 1000.0d0
    D_o_m = D_outer_mm / 1000.0d0
    L_e_m = Le_mm / 1000.0d0
    L_a_m = La_mm / 1000.0d0
    L_c_m = Lc_mm / 1000.0d0

    L_eff_m = L_a_m + (L_e_m + L_c_m) / 2.0d0
    A_v_m2 = PI * (D_i_m**2) / 4.0d0
    A_e_m2 = PI * D_i_m * L_e_m
    A_c_m2 = PI * D_i_m * L_c_m

    ! 1. Entrainment / Flooding Limit (Wallis / Faghri correlation)
    Q_flood_W = 1.8d0 * A_v_m2 * (hfg_kJkg * 1000.0d0) * ((G_ACC * sigma_Nm * (rho_L - rho_v))**0.25d0) * &
                ((rho_v**(-0.25d0) + rho_L**(-0.25d0))**(-2.0d0))

    ! 2. Boiling Dryout Limit
    Q_boiling_W = 0.18d0 * A_e_m2 * (hfg_kJkg * 1000.0d0) * sqrt(rho_v) * ((G_ACC * sigma_Nm * (rho_L - rho_v))**0.25d0)

    ! 3. Sonic Limit
    Q_sonic_W = 0.474d0 * A_v_m2 * (hfg_kJkg * 1000.0d0) * sqrt(1.33d0 * rho_v * (Pv_bar * 1.0d5))

    Q_max_allow_W = min(Q_flood_W, min(Q_boiling_W, Q_sonic_W))

    if (Q_max_allow_W == Q_flood_W) then
        limit_type = 'FLOODING / ENTRAINMENT LIMIT'
    else if (Q_max_allow_W == Q_boiling_W) then
        limit_type = 'BOILING BURNOUT LIMIT'
    else
        limit_type = 'SONIC CHOKING LIMIT'
    end if

    ! Overall thermal drop ΔT approx 2 - 8 K across thermosiphon
    delta_T_pipe_C = (Q_heat_W / (2500.0d0 * A_e_m2)) + (Q_heat_W / (4000.0d0 * A_c_m2)) + 0.5d0

    ! Effective thermal conductivity k_eff [W/(m.K)]
    k_eff_WmK = (Q_heat_W * L_eff_m) / (A_v_m2 * max(0.1d0, delta_T_pipe_C))

    if (Q_heat_W > Q_max_allow_W) then
        status_str = 'WARNING: EXCEEDS LIMIT (BURNOUT RISK)'
    else if (Q_heat_W > 0.70d0 * Q_max_allow_W) then
        status_str = 'SAFE (APPROACHING PEAK CAPACITY)'
    else
        status_str = 'OPTIMAL SAFE OPERATION'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — TWO-PHASE CLOSED THERMOSIPHON / HEAT PIPE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Vapor Sat Temp / Pressure    = ', Tv_sat_C, ' deg C / ', Pv_bar, ' bar abs'
    write(*,'(A,F10.1,A,F10.1,A)') 'Pipe OD / Effective Length   = ', D_outer_mm, ' mm / ', L_eff_m*1000.0d0, ' mm'
    write(*,'(A,A)')               'Dominant Critical Limit      = ', trim(limit_type)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.1,A)')  'HEAT LOAD CARRIED (Q)        = ', Q_heat_W, ' W'
    write(*,'(A,F10.1,A)')  'MAX ALLOWABLE CAPACITY (Qmax)= ', Q_max_allow_W, ' W'
    write(*,'(A,F10.1,A)')  'Flooding / Entrainment Limit = ', Q_flood_W, ' W'
    write(*,'(A,F10.1,A)')  'Boiling Pool Burnout Limit   = ', Q_boiling_W, ' W'
    write(*,'(A,F10.2,A)')  'End-to-End Temp Drop (ΔT)    = ', delta_T_pipe_C, ' deg C'
    write(*,'(A,F10.0,A)')  'EFFECTIVE THERMAL CONDUCTIVITY= ', k_eff_WmK, ' W/(m.K) (Superconductor)'
    write(*,'(A,A)')        'OPERATIONAL SAFETY STATUS    = ', trim(status_str)
    write(*,'(A)') '============================================================'

end program two_phase_closed_thermosiphon_heat_pipe

💻 How to Compile & Run Locally

1. Compilation (GNU Fortran / Intel oneAPI):

gfortran -O3 two_phase_closed_thermosiphon_heat_pipe.f90 -o two_phase_closed_thermosiphon_heat_pipe

2. Execution with input.txt redirection:

two_phase_closed_thermosiphon_heat_pipe < input.txt

📄 Sample input.txt File Structure

Sample Data:

Parameter Description:

Fluid (1=Water, 2=Ethanol, 3=Methanol, 4=Acetone, 5=Ammonia)\nPipe Material (1=Copper, 2=Aluminum, 3=Stainless)\nHeat Load [W]\nVapor Saturation Temp [°C]\nPipe OD [mm]\nPipe ID [mm]\nEvaporator Length [mm]\nAdiabatic Length [mm]\nCondenser Length [mm]\nFill Ratio [%]