💻 Fortran Source Code Library

! =========================================================================
! Source File: heat_sink_calculator.f90
! =========================================================================

program Heat_Sink_Calculator
    implicit none

    ! Inputs
    integer :: type_id          ! 1=Plate fin, 2=Pin fin array
    real(8) :: width            ! mm (base plate width W)
    real(8) :: depth            ! mm (base plate depth D)
    real(8) :: t_base           ! mm (base plate thickness)
    integer :: n_fins           ! total fins/pins
    real(8) :: h_fin            ! mm (fin height H)
    real(8) :: t_fin            ! mm (fin thickness t or pin diameter d)
    real(8) :: s_fin            ! mm (fin spacing, read but calculated dynamically)
    integer :: mat_id           ! 1=Al, 2=Cu, 3=Steel, 4=Brass, 5=Cast Iron, 6=Custom
    real(8) :: k_custom         ! W/mK
    integer :: flow_mode        ! 1=Forced Convection, 2=Natural Convection
    real(8) :: v_in            ! m/s (inlet velocity)
    real(8) :: power            ! W (heat Q)
    real(8) :: t_amb            ! degC (ambient temp)
    integer :: h_custom_flag    ! 0=Auto-compute, 1=Manual override
    real(8) :: h_manual         ! W/m2K (override value)

    ! SI Variables (using unique case-insensitive prefixes)
    real(8) :: w_base_m, d_depth_m, t_plate_m, h_fin_m, t_fin_m, d_pin_m
    real(8) :: k_mat, h_coeff
    real(8) :: T_base_temp, T_base_old_temp, T_film_temp, T_film_K_temp
    real(8) :: rho, mu, nu, Cp, k_air, Pr, alpha, beta
    real(8) :: g, PI
    integer :: iter, max_iter

    ! Geometry & Performance variables
    real(8) :: b_spacing, U_ch, Re_b, Re_b_star, Nu_b, term1, term2
    integer :: N_x, N_y
    real(8) :: S_T, S_L, S_pitch_x, V_max, Re_d_max, Nu_d, Ra_d, Ra_S, S_mean
    real(8) :: m_val, H_c, eta_f, A_f, A_b, A_tot, eta_o, R_base, R_hs
    real(8) :: alpha_star, f_D, f_lam, f_turb, sigma, K_c, K_e, D_h, dp_val
    real(8) :: Ra_H, S_opt, Re_D_forced, Ra_b_star

    ! Parameters
    g = 9.81d0
    PI = 3.14159265358979323846d0

    ! ---------------------------------------------------------
    ! 1. READ INPUT PARAMETERS
    ! ---------------------------------------------------------
    read(*,*) type_id
    read(*,*) width
    read(*,*) depth
    read(*,*) t_base
    read(*,*) n_fins
    read(*,*) h_fin
    read(*,*) t_fin
    read(*,*) s_fin
    read(*,*) mat_id
    read(*,*) k_custom
    read(*,*) flow_mode
    read(*,*) v_in
    read(*,*) power
    read(*,*) t_amb
    read(*,*) h_custom_flag
    read(*,*) h_manual

    ! Convert to SI
    w_base_m = width / 1000.0d0
    d_depth_m = depth / 1000.0d0
    t_plate_m = t_base / 1000.0d0
    h_fin_m = h_fin / 1000.0d0
    t_fin_m = t_fin / 1000.0d0
    d_pin_m = t_fin / 1000.0d0

    ! Material presets
    select case(mat_id)
    case(1)
        k_mat = 200.0d0  ! Aluminum (6061-T6)
    case(2)
        k_mat = 385.0d0  ! Copper (pure)
    case(3)
        k_mat = 50.0d0   ! Steel (low-carbon)
    case(4)
        k_mat = 110.0d0  ! Brass
    case(5)
        k_mat = 52.0d0   ! Cast Iron
    case(6)
        k_mat = k_custom
    case default
        k_mat = 200.0d0
    end select

    ! Initial guess for base temperature
    T_base_temp = t_amb + 20.0d0
    max_iter = 10

    ! ---------------------------------------------------------
    ! 2. ITERATIVE LOOP TO CONVERGE T_BASE AND FLUID PROPERTIES
    ! ---------------------------------------------------------
    do iter = 1, max_iter
        T_base_old_temp = T_base_temp
        T_film_temp = (T_base_temp + t_amb) / 2.0d0
        T_film_K_temp = T_film_temp + 273.15d0

        ! Air properties at film temperature
        rho = 101325.0d0 / (287.05d0 * T_film_K_temp)
        mu = 1.458d-6 * (T_film_K_temp**1.5d0) / (T_film_K_temp + 110.4d0)
        nu = mu / rho
        Cp = 1005.0d0
        k_air = 0.0263d0 * (T_film_K_temp / 293.15d0)**0.85d0
        Pr = 0.71d0
        alpha = k_air / (rho * Cp)
        beta = 1.0d0 / T_film_K_temp

        ! Convection coefficient calculation
        if (h_custom_flag == 1) then
            h_coeff = h_manual
        else
            if (type_id == 1) then
                ! Plate Fin
                b_spacing = (w_base_m - n_fins * t_fin_m) / max(real(n_fins - 1, 8), 1.0d0)
                if (b_spacing < 1.0d-4) b_spacing = 1.0d-4

                if (flow_mode == 1) then
                    ! Forced Convection (Teertstra correlation)
                    U_ch = v_in * w_base_m / max((n_fins - 1) * b_spacing, 1.0d-6)
                    Re_b = U_ch * b_spacing / nu
                    Re_b_star = Re_b * b_spacing / d_depth_m
                    
                    term1 = (Re_b_star * Pr / 2.0d0)**(-3.0d0)
                    term2 = (0.664d0 * sqrt(Re_b_star) * (Pr**(1.0d0/3.0d0)) * &
                             sqrt(1.0d0 + 3.65d0 / sqrt(max(Re_b_star, 1.0d-8))))**(-3.0d0)
                    Nu_b = (term1 + term2)**(-1.0d0/3.0d0)
                    h_coeff = Nu_b * k_air / b_spacing
                else
                    ! Natural Convection (Elenbaas correlation)
                    Ra_b_star = (g * beta * abs(T_base_temp - t_amb) * (b_spacing**3.0d0) / &
                                 (nu * alpha)) * (b_spacing / h_fin_m)
                    
                    if (Ra_b_star > 1.0d-5) then
                        Nu_b = ((576.0d0 / (Ra_b_star**2)) + (2.873d0 / sqrt(Ra_b_star)))**(-0.5d0)
                    else
                        Nu_b = Ra_b_star / 24.0d0
                    end if
                    h_coeff = Nu_b * k_air / b_spacing
                end if
            else
                ! Pin Fin Array
                ! Calculate Nx, Ny
                N_x = nint(sqrt(real(n_fins, 8) * w_base_m / d_depth_m))
                if (N_x < 2) N_x = 2
                N_y = nint(real(n_fins, 8) / real(N_x, 8))
                if (N_y < 1) N_y = 1
                
                S_T = (w_base_m - N_x * d_pin_m) / max(real(N_x - 1, 8), 1.0d0)
                S_L = (d_depth_m - N_y * d_pin_m) / max(real(N_y - 1, 8), 1.0d0)
                if (S_T < 1.0d-4) S_T = 1.0d-4
                if (S_L < 1.0d-4) S_L = 1.0d-4

                if (flow_mode == 1) then
                    ! Forced Convection (Zukauskas correlation)
                    S_pitch_x = S_T + d_pin_m
                    V_max = v_in * S_pitch_x / max(S_T, 1.0d-6)
                    Re_d_max = V_max * d_pin_m / nu
                    
                    if (Re_d_max < 100.0d0) then
                        Nu_d = 0.90d0 * (Re_d_max**0.4d0) * (Pr**0.36d0)
                    else if (Re_d_max < 1000.0d0) then
                        Nu_d = 0.51d0 * (Re_d_max**0.5d0) * (Pr**0.36d0)
                    else if (Re_d_max < 200000.0d0) then
                        Nu_d = 0.35d0 * (Re_d_max**0.6d0) * (Pr**0.36d0)
                    else
                        Nu_d = 0.021d0 * (Re_d_max**0.84d0) * (Pr**0.36d0)
                    end if
                    h_coeff = Nu_d * k_air / d_pin_m
                else
                    ! Natural Convection (Vertical cylinder correlation)
                    Ra_d = g * beta * abs(T_base_temp - t_amb) * (d_pin_m**3.0d0) / (nu * alpha)
                    if (Ra_d > 0.0d0) then
                        if (Ra_d < 1.0d9) then
                            Nu_d = 0.54d0 * (Ra_d**0.25d0)
                        else
                            Nu_d = 0.13d0 * (Ra_d**(1.0d0/3.0d0))
                        end if
                    else
                        Nu_d = 0.5d0
                    end if
                    h_coeff = Nu_d * k_air / d_pin_m
                end if
            end if
        end if

        ! Ensure h is non-zero
        if (h_coeff < 0.1d0) h_coeff = 0.1d0

        ! ---------------------------------------------------------
        ! 3. FIN EFFICIENCY AND THERMAL RESISTANCE
        ! ---------------------------------------------------------
        if (type_id == 1) then
            ! Plate Fin efficiency
            m_val = sqrt(2.0d0 * h_coeff / max(k_mat * t_fin_m, 1.0d-8))
            H_c = h_fin_m + t_fin_m / 2.0d0
            eta_f = tanh(m_val * H_c) / (m_val * H_c)
            A_f = 2.0d0 * H_c * d_depth_m
            A_b = (w_base_m - n_fins * t_fin_m) * d_depth_m
            if (A_b < 0.0d0) A_b = 0.0d0
            A_tot = n_fins * A_f + A_b
            eta_o = 1.0d0 - (real(n_fins, 8) * A_f / A_tot) * (1.0d0 - eta_f)
        else
            ! Pin Fin array efficiency
            N_x = nint(sqrt(real(n_fins, 8) * w_base_m / d_depth_m))
            if (N_x < 2) N_x = 2
            N_y = nint(real(n_fins, 8) / real(N_x, 8))
            if (N_y < 1) N_y = 1
            
            m_val = sqrt(4.0d0 * h_coeff / max(k_mat * d_pin_m, 1.0d-8))
            H_c = h_fin_m + d_pin_m / 4.0d0
            eta_f = tanh(m_val * H_c) / (m_val * H_c)
            A_f = PI * d_pin_m * H_c
            A_b = w_base_m * d_depth_m - real(N_x * N_y, 8) * (PI * d_pin_m**2 / 4.0d0)
            if (A_b < 0.0d0) A_b = 0.0d0
            A_tot = real(N_x * N_y, 8) * A_f + A_b
            eta_o = 1.0d0 - (real(N_x * N_y, 8) * A_f / A_tot) * (1.0d0 - eta_f)
        end if

        R_base = t_plate_m / max(k_mat * w_base_m * d_depth_m, 1.0d-8)
        R_hs = R_base + 1.0d0 / max(eta_o * h_coeff * A_tot, 1.0d-12)
        T_base_temp = t_amb + power * R_hs

        ! Convergence check
        if (abs(T_base_temp - T_base_old_temp) < 1.0d-4) exit
    end do

    ! ---------------------------------------------------------
    ! 4. PRESSURE DROP COMPUTATION
    ! ---------------------------------------------------------
    dp_val = 0.0d0
    if (flow_mode == 1) then
        if (type_id == 1) then
            ! Plate Fin channel flow pressure drop
            b_spacing = (w_base_m - n_fins * t_fin_m) / max(real(n_fins - 1, 8), 1.0d0)
            if (b_spacing < 1.0d-4) b_spacing = 1.0d-4

            U_ch = v_in * w_base_m / max((n_fins - 1) * b_spacing, 1.0d-6)
            D_h = 2.0d0 * b_spacing * h_fin_m / (b_spacing + h_fin_m)
            Re_b = U_ch * D_h / nu

            ! Aspect ratio
            alpha_star = b_spacing / h_fin_m
            if (alpha_star > 1.0d0) alpha_star = 1.0d0 / alpha_star

            ! Friction factor
            if (Re_b < 2300.0d0) then
                f_lam = 96.0d0 * (1.0d0 - 1.3553d0 * alpha_star + 1.9467d0 * (alpha_star**2) - &
                                  1.7012d0 * (alpha_star**3) + 0.9564d0 * (alpha_star**4) - &
                                  0.2537d0 * (alpha_star**5)) / max(Re_b, 1.0d0)
                f_D = f_lam
            else
                ! Blasius turbulent friction factor
                f_turb = 0.316d0 / (max(Re_b, 2300.0d0)**0.25d0)
                f_D = f_turb
            end if

            sigma = b_spacing / (b_spacing + t_fin_m)
            K_c = 0.8d0 - 0.4d0 * (sigma**2)
            K_e = 1.2d0 - 1.6d0 * sigma + 0.4d0 * (sigma**2)
            dp_val = (K_c + K_e + f_D * d_depth_m / D_h) * 0.5d0 * rho * (U_ch**2)
        else
            ! Pin Fin array cross-flow pressure drop
            N_x = nint(sqrt(real(n_fins, 8) * w_base_m / d_depth_m))
            if (N_x < 2) N_x = 2
            N_y = nint(real(n_fins, 8) / real(N_x, 8))
            if (N_y < 1) N_y = 1
            
            S_T = (w_base_m - N_x * d_pin_m) / max(real(N_x - 1, 8), 1.0d0)
            if (S_T < 1.0d-4) S_T = 1.0d-4

            S_pitch_x = S_T + d_pin_m
            V_max = v_in * S_pitch_x / max(S_T, 1.0d-6)
            Re_d_max = V_max * d_pin_m / nu

            f_D = 0.2d0 + 50.0d0 / max(Re_d_max, 1.0d-2) + 5.0d0 / max(sqrt(Re_d_max), 1.0d-2)
            dp_val = real(N_y, 8) * f_D * 0.5d0 * rho * (V_max**2)
        end if
    end if

    ! ---------------------------------------------------------
    ! 5. OPTIMAL FIN SPACING
    ! ---------------------------------------------------------
    S_opt = 0.0d0
    if (type_id == 1) then
        if (flow_mode == 2) then
            ! Natural Convection optimal spacing (Bar-Cohen & Rohsenow)
            Ra_H = g * beta * abs(T_base_temp - t_amb) * (h_fin_m**3.0d0) / (nu * alpha)
            if (Ra_H > 1.0d-2) then
                S_opt = 2.714d0 * h_fin_m / (Ra_H**0.25d0) * 1000.0d0 ! mm
            else
                S_opt = 0.0d0
            end if
        else
            ! Forced Convection optimal spacing (Bejan parallel plates)
            Re_D_forced = v_in * d_depth_m / nu
            if (Re_D_forced > 1.0d-2) then
                S_opt = 2.73d0 * d_depth_m / (sqrt(Re_D_forced) * (Pr**0.25d0)) * 1000.0d0 ! mm
            else
                S_opt = 0.0d0
            end if
        end if
    end if

    ! ---------------------------------------------------------
    ! 6. OUTPUT RESULTS REPORT
    ! ---------------------------------------------------------
    write(*,*) '================================================================'
    write(*,*) '                    INPUT PARAMETERS'
    write(*,*) '================================================================'
    write(*,*)
    write(*,'(A,I2)')      ' Heat Sink Type ID:            ', type_id
    write(*,'(A,F10.2,A)') ' Base Plate Width (W):         ', width, ' mm'
    write(*,'(A,F10.2,A)') ' Base Plate Depth (D):         ', depth, ' mm'
    write(*,'(A,F10.2,A)') ' Base Plate Thickness (t):     ', t_base, ' mm'
    write(*,'(A,I5)')      ' Number of Fins/Pins (N):      ', n_fins
    write(*,'(A,F10.2,A)') ' Fin Height (H):               ', h_fin, ' mm'
    write(*,'(A,F10.2,A)') ' Fin Thickness/Diameter:       ', t_fin, ' mm'
    write(*,'(A,I2)')      ' Flow Convection Mode:         ', flow_mode
    if (flow_mode == 1) then
        write(*,'(A,F10.2,A)') ' Inlet Air Velocity (V_in):    ', v_in, ' m/s'
    end if
    write(*,'(A,F10.2,A)') ' Heat Dissipation load (Q):    ', power, ' W'
    write(*,'(A,F10.2,A)') ' Ambient Temperature (T_inf):  ', t_amb, ' degC'
    write(*,*)

    write(*,*) '================================================================'
    write(*,*) '                   THERMAL HYDRAULIC RESULTS'
    write(*,*) '================================================================'
    write(*,*)
    write(*,'(A,F12.4,A)') ' Base Temperature (T_base):    ', T_base_temp, ' degC'
    write(*,'(A,F12.6,A)') ' Thermal Resistance (R_hs):    ', R_hs, ' degC/W'
    write(*,'(A,F12.4,A)') ' Convection Coeff. (h):        ', h_coeff, ' W/m2.K'
    write(*,'(A,F12.2,A)') ' Overall Surface Area (A_tot): ', A_tot * 10000.0d0, ' cm2'
    write(*,'(A,F12.2,A)') ' Fin Efficiency (eta_f):       ', eta_f * 100.0d0, ' %'
    write(*,'(A,F12.2,A)') ' Overall efficiency (eta_o):   ', eta_o * 100.0d0, ' %'
    if (flow_mode == 1) then
        write(*,'(A,F12.4,A)') ' Pressure Drop (Delta P):      ', dp_val, ' Pa'
    else
        write(*,'(A,F12.4,A)') ' Pressure Drop (Delta P):      ', 0.0d0, ' Pa (Natural)'
    end if
    if (type_id == 1) then
        write(*,'(A,F12.3,A)') ' Calculated Fin Spacing (S):   ', b_spacing * 1000.0d0, ' mm'
        if (S_opt > 0.0d0) then
            write(*,'(A,F12.3,A)') ' Optimal Fin Spacing (S_opt):  ', S_opt, ' mm'
        end if
    else
        write(*,'(A,F12.3,A)') ' Transverse spacing (S_T):     ', S_T * 1000.0d0, ' mm'
        write(*,'(A,F12.3,A)') ' Longitudinal spacing (S_L):   ', S_L * 1000.0d0, ' mm'
    end if
    write(*,*)

    write(*,*) '================================================================'
    write(*,*) '                MATERIAL PERFORMANCE COMPARISON'
    write(*,*) '================================================================'
    write(*,*) '  (Same geometry, flow mode, power, and ambient temperature)'
    write(*,*)
    call compare_mat_hs('Aluminum (6061-T6)', k_mat, R_base, A_tot, eta_f, H_c, h_coeff)
    call compare_mat_hs('Copper (pure)     ', k_mat, R_base, A_tot, eta_f, H_c, h_coeff)
    call compare_mat_hs('Steel (low-carbon)', k_mat, R_base, A_tot, eta_f, H_c, h_coeff)
    call compare_mat_hs('Brass (Cu65/Zn35) ', k_mat, R_base, A_tot, eta_f, H_c, h_coeff)
    call compare_mat_hs('Cast Iron         ', k_mat, R_base, A_tot, eta_f, H_c, h_coeff)
    write(*,*)
    write(*,*) '================================================================'
    write(*,*) '               CALCULATION COMPLETE'
    write(*,*) '================================================================'

contains

    subroutine compare_mat_hs(name, ref_k, ref_R_base, ref_A_tot, ref_eta_f, ref_H_c, h_cf)
        implicit none
        character(len=*), intent(in) :: name
        real(8), intent(in) :: ref_k, ref_R_base, ref_A_tot, ref_eta_f, ref_H_c, h_cf
        real(8) :: k_m, m_m, eta_fm, eta_om, R_basem, R_hsm, T_basem

        ! Lookup thermal conductivity
        select case(name)
        case('Aluminum (6061-T6)')
            k_m = 200.0d0
        case('Copper (pure)     ')
            k_m = 385.0d0
        case('Steel (low-carbon)')
            k_m = 50.0d0
        case('Brass (Cu65/Zn35) ')
            k_m = 110.0d0
        case('Cast Iron         ')
            k_m = 52.0d0
        case default
            k_m = 200.0d0
        end select

        ! Recalculate efficiency for this material
        if (type_id == 1) then
            m_m = sqrt(2.0d0 * h_cf / max(k_m * t_fin_m, 1.0d-8))
            eta_fm = tanh(m_m * ref_H_c) / (m_m * ref_H_c)
            eta_om = 1.0d0 - (real(n_fins, 8) * A_f / ref_A_tot) * (1.0d0 - eta_fm)
        else
            m_m = sqrt(4.0d0 * h_cf / max(k_m * d_pin_m, 1.0d-8))
            eta_fm = tanh(m_m * ref_H_c) / (m_m * ref_H_c)
            eta_om = 1.0d0 - (real(N_x * N_y, 8) * A_f / ref_A_tot) * (1.0d0 - eta_fm)
        end if

        R_basem = t_plate_m / max(k_m * w_base_m * d_depth_m, 1.0d-8)
        R_hsm = R_basem + 1.0d0 / max(eta_om * h_cf * ref_A_tot, 1.0d-12)
        T_basem = t_amb + power * R_hsm

        write(*,'(3X,A20,A,F7.4,A,F7.2,A)') name, ': R_hs = ', R_hsm, &
            ' degC/W  →  T_base = ', T_basem, ' degC'
    end subroutine compare_mat_hs

end program Heat_Sink_Calculator