Cavity Enclosure Natural Convection (Catton)

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

Model natural convection in differentially heated vertical cavity enclosures (double glazing, solar gaps), computing Rayleigh number, average Nusselt, total heat flow Q, and boundary layer thickness.

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

📁 calcul/Convection / natural_convection_cavity_enclosure.f90

program natural_convection_cavity_enclosure
    implicit none
    integer :: iostat_val
    double precision :: H_height_m, W_width_m, L_depth_m, Th_hot_C, Tc_cold_C
    double precision :: rho_dens, mu_cP, k_cond, cp_spec, beta_exp
    double precision :: mu_Pas, nu_visc, alpha_diff, DeltaT, Ra_W, Pr_num, AspectRatio
    double precision :: Nu_avg, h_avg, Q_dot_W, delta_bl_mm, u_buoy_ms
    double precision, parameter :: G_GRAV = 9.80665d0

    ! Read inputs
    read(*,*,iostat=iostat_val) H_height_m    ! Cavity Height [m] (e.g. 1.0)
    read(*,*,iostat=iostat_val) W_width_m     ! Cavity Spacing Width [m] (e.g. 0.05)
    read(*,*,iostat=iostat_val) L_depth_m     ! Cavity Depth [m] (e.g. 1.0)
    read(*,*,iostat=iostat_val) Th_hot_C      ! Hot Wall Temp [deg C] (e.g. 50.0)
    read(*,*,iostat=iostat_val) Tc_cold_C     ! Cold Wall Temp [deg C] (e.g. 10.0)
    read(*,*,iostat=iostat_val) rho_dens      ! Fluid Density [kg/m3] (e.g. 1.164 for air)
    read(*,*,iostat=iostat_val) mu_cP         ! Fluid Viscosity [cP] (e.g. 0.0185 for air)
    read(*,*,iostat=iostat_val) k_cond        ! Thermal Conductivity [W/(m.K)] (e.g. 0.026)
    read(*,*,iostat=iostat_val) cp_spec       ! Specific Heat [J/(kg.K)] (e.g. 1005.0)
    read(*,*,iostat=iostat_val) beta_exp      ! Volumetric Expansion Coeff [1/K] (e.g. 0.00333)

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

    if (H_height_m <= 0.0d0 .or. W_width_m <= 0.0d0 .or. L_depth_m <= 0.0d0) then
        write(*,*) 'ERROR: Geometry dimensions must be positive.'
        stop
    end if

    DeltaT = abs(Th_hot_C - Tc_cold_C)
    if (DeltaT <= 0.0d0) DeltaT = 1.0d-3

    mu_Pas = mu_cP * 1.0d-3
    nu_visc = mu_Pas / rho_dens
    alpha_diff = k_cond / (rho_dens * cp_spec)
    Pr_num = nu_visc / alpha_diff
    AspectRatio = H_height_m / W_width_m

    Ra_W = (G_GRAV * beta_exp * DeltaT * (W_width_m**3)) / (nu_visc * alpha_diff)

    ! Catton / Berkovsky-Polevikov Correlation for Vertical Enclosure
    if (Ra_W < 1000.0d0) then
        Nu_avg = 1.0d0 ! Pure conduction
    else if (Ra_W < 1.0d4) then
        Nu_avg = 0.18d0 * ((Ra_W * Pr_num / (0.2d0 + Pr_num))**0.29d0)
    else
        Nu_avg = 0.22d0 * ((Ra_W * Pr_num / (0.2d0 + Pr_num))**0.28d0) * (AspectRatio**(-0.25d0))
    end if
    if (Nu_avg < 1.0d0) Nu_avg = 1.0d0

    h_avg = (Nu_avg * k_cond) / W_width_m
    Q_dot_W = h_avg * (H_height_m * L_depth_m) * DeltaT
    delta_bl_mm = (W_width_m / (max(1.0d0, Ra_W)**0.25d0)) * 1000.0d0
    u_buoy_ms = sqrt(G_GRAV * beta_exp * DeltaT * H_height_m)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — CAVITY NATURAL CONVECTION ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Cavity H x W (Aspect A)   = ', H_height_m, ' m x ', W_width_m, ' m (A = ', AspectRatio, ')'
    write(*,'(A,F10.1,A,F10.1,A)') 'Hot Wall Th / Cold Wall Tc= ', Th_hot_C, ' deg C / ', Tc_cold_C, ' deg C'
    write(*,'(A,E12.4,A,F10.2)')   'Cavity Rayleigh Ra_W / Pr = ', Ra_W, ' / ', Pr_num
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2)')           'AVERAGE NUSSELT (Nu_avg)  = ', Nu_avg
    write(*,'(A,F10.2,A)')  'Heat Transfer Coeff h_avg = ', h_avg, ' W/(m2.K)'
    write(*,'(A,F10.2,A,F10.2,A)') 'TOTAL HEAT FLOW (Q_dot)   = ', Q_dot_W, ' W (', Q_dot_W/1000.0d0, ' kW)'
    write(*,'(A,F10.2,A,F10.3,A)') 'BL Thickness / Max Buoy U = ', delta_bl_mm, ' mm / ', u_buoy_ms, ' m/s'
    write(*,'(A)') '============================================================'

end program natural_convection_cavity_enclosure

💻 How to Compile & Run Locally

1. Compilation (GNU Fortran / Intel oneAPI):

gfortran -O3 natural_convection_cavity_enclosure.f90 -o natural_convection_cavity_enclosure

2. Execution with input.txt redirection:

natural_convection_cavity_enclosure < input.txt

📄 Sample input.txt File Structure

Sample Data:

Parameter Description:

Cavity Height H [m]\nCavity Spacing W [m]\nCavity Depth L [m]\nHot Wall Temp [°C]\nCold Wall Temp [°C]\nDensity [kg/m³]\nViscosity [cP]\nConductivity [W/(m·K)]\nSpecific Heat [J/(kg·K)]\nExpansion Coeff beta [1/K]