💻 Fortran Source Code Library

! =========================================================================
! Source File: pin_fin_calculator.f90
! =========================================================================

program Pin_Fin_Spine_Calculator
    implicit none

    ! Inputs
    integer :: shape_id         ! 1=Cylindrical, 2=Conical, 3=Parabolic, 4=Rectangular
    real(8) :: d_base           ! mm (diameter for circular, width for rectangular)
    real(8) :: t_base_rect      ! mm (thickness for rectangular, not used for others)
    real(8) :: length           ! mm
    real(8) :: k_cond           ! W/m.K
    real(8) :: h_conv           ! W/m2.K
    real(8) :: T_base_temp      ! degC
    real(8) :: T_ambient        ! degC
    integer :: tip_cond         ! 1=Adiabatic, 2=Convective, 3=Prescribed T, 4=Infinite
    real(8) :: T_tip            ! degC (only if tip_cond = 3)
    real(8) :: h_tip            ! W/m2.K (only if tip_cond = 2)

    ! Internal variables
    real(8) :: d_m, t_m, L_m
    real(8) :: Ac, P, As, Atot
    real(8) :: m_val, mL, theta_b, theta_L
    real(8) :: Q_fin, eta_f, epsilon_f, Biot_f
    real(8) :: Qmax
    real(8) :: theta_x, T_x, x_val, z_val
    real(8), parameter :: PI = 3.14159265358979323846d0
    integer :: i, n_points
    real(8) :: I0_val, I1_val, I0_base, I1_base, I2_base
    real(8) :: a_param, opt_L
    real(8) :: arg_base, arg_z

    ! ---------------------------------------------------------
    ! 1. READ INPUT PARAMETERS
    ! ---------------------------------------------------------
    read(*,*) shape_id
    read(*,*) d_base
    read(*,*) t_base_rect
    read(*,*) length
    read(*,*) k_cond
    read(*,*) h_conv
    read(*,*) T_base_temp
    read(*,*) T_ambient
    read(*,*) tip_cond
    read(*,*) T_tip
    read(*,*) h_tip

    ! Conversions to SI
    d_m = d_base / 1000.0d0
    t_m = t_base_rect / 1000.0d0
    L_m = length / 1000.0d0
    theta_b = T_base_temp - T_ambient

    ! ---------------------------------------------------------
    ! 2. GEOMETRY AND PARAMETER EVALUATION
    ! ---------------------------------------------------------
    select case(shape_id)
    case(1) ! Cylindrical Spine
        Ac = PI * d_m**2 / 4.0d0
        P = PI * d_m
        As = P * L_m
        if (tip_cond == 2) then
            Atot = As + Ac
        else
            Atot = As
        end if
        m_val = sqrt(4.0d0 * h_conv / (k_cond * d_m))
    case(4) ! Rectangular Spine
        Ac = d_m * t_m
        P = 2.0d0 * (d_m + t_m)
        As = P * L_m
        if (tip_cond == 2) then
            Atot = As + Ac
        else
            Atot = As
        end if
        m_val = sqrt((h_conv * P) / (k_cond * Ac))
    case(2) ! Conical Spine
        Ac = PI * d_m**2 / 4.0d0
        P = PI * d_m
        As = PI * (d_m / 2.0d0) * sqrt(L_m**2 + (d_m / 2.0d0)**2)
        Atot = As
        m_val = sqrt(4.0d0 * h_conv / (k_cond * d_m))
    case(3) ! Parabolic Spine (Convex)
        Ac = PI * d_m**2 / 4.0d0
        P = PI * d_m
        a_param = d_m / L_m
        if (a_param > 1d-6) then
            As = PI * (d_m / 2.0d0) * L_m * ( ((2.0d0 * a_param**2 + 1.0d0) * &
                 sqrt(a_param**2 + 1.0d0)) / (8.0d0 * a_param**2) - &
                 log(a_param + sqrt(a_param**2 + 1.0d0)) / (8.0d0 * a_param**3) )
        else
            As = PI * d_m * L_m / 2.0d0
        end if
        Atot = As
        m_val = sqrt(4.0d0 * h_conv / (k_cond * d_m))
    end select

    mL = m_val * L_m
    Qmax = h_conv * Atot * theta_b

    ! ---------------------------------------------------------
    ! 3. EFFICIENCY AND HEAT TRANSFER DISSIPATION
    ! ---------------------------------------------------------
    if (shape_id == 1 .or. shape_id == 4) then
        ! Uniform Spine
        select case(tip_cond)
        case(1) ! Adiabatic
            eta_f = tanh(mL) / mL
            Q_fin = eta_f * Qmax
        case(2) ! Convective
            Q_fin = sqrt(h_conv * P * k_cond * Ac) * theta_b * &
                    (sinh(mL) + (h_tip / (m_val * k_cond)) * cosh(mL)) / &
                    (cosh(mL) + (h_tip / (m_val * k_cond)) * sinh(mL))
            eta_f = Q_fin / (h_conv * As * theta_b + h_tip * Ac * theta_b)
        case(3) ! Prescribed T
            theta_L = T_tip - T_ambient
            Q_fin = sqrt(h_conv * P * k_cond * Ac) * theta_b * &
                    (cosh(mL) - theta_L / theta_b) / sinh(mL)
            eta_f = Q_fin / Qmax
        case(4) ! Infinite
            Q_fin = sqrt(h_conv * P * k_cond * Ac) * theta_b
            eta_f = 1.0d0 / mL
        end select
    else if (shape_id == 2) then
        ! Conical Spine
        arg_base = 2.0d0 * mL
        call bessel_I1(arg_base, I1_base)
        call bessel_I2(arg_base, I2_base)
        if (I1_base > 1d-30) then
            eta_f = (2.0d0 / mL) * (I2_base / I1_base)
        else
            eta_f = 1.0d0
        end if
        Q_fin = eta_f * h_conv * As * theta_b
    else if (shape_id == 3) then
        ! Parabolic Spine (Convex)
        arg_base = (4.0d0 / 3.0d0) * mL
        call bessel_I0(arg_base, I0_base)
        call bessel_I1(arg_base, I1_base)
        if (I0_base > 1d-30) then
            eta_f = (1.5d0 / mL) * (I1_base / I0_base)
        else
            eta_f = 1.0d0
        end if
        Q_fin = eta_f * h_conv * As * theta_b
    end if

    ! Limit efficiency to physical bounds
    if (eta_f > 1.0d0) eta_f = 1.0d0
    if (eta_f < 0.0d0) eta_f = 0.0d0

    Biot_f = h_conv * (Ac / P) / k_cond
    epsilon_f = Q_fin / max(h_conv * Ac * theta_b, 1.0d-30)

    ! Write Output Summary
    write(*,*) '================================================================'
    write(*,*) '                    INPUT PARAMETERS'
    write(*,*) '================================================================'
    write(*,*)
    write(*,'(A,I2)')      ' Shape ID:                     ', shape_id
    write(*,'(A,I2)')      ' Tip Condition ID:             ', tip_cond
    write(*,'(A,F10.4,A)') ' Base diameter/width (d/w):    ', d_base, ' mm'
    if (shape_id == 4) then
        write(*,'(A,F10.4,A)') ' Thickness/height (t):         ', t_base_rect, ' mm'
    end if
    write(*,'(A,F10.4,A)') ' Fin Length (L):               ', length, ' mm'
    write(*,'(A,F10.3,A)') ' Thermal Conductivity (k):     ', k_cond, ' W/m.K'
    write(*,'(A,F10.3,A)') ' Convection Coefficient (h):   ', h_conv, ' W/m2.K'
    write(*,'(A,F10.2,A)') ' Base Temperature (T_b):       ', T_base_temp, ' deg-C'
    write(*,'(A,F10.2,A)') ' Ambient Temperature (T_inf):  ', T_ambient, ' deg-C'
    write(*,*)

    write(*,*) '================================================================'
    write(*,*) '                   FIN PARAMETERS'
    write(*,*) '================================================================'
    write(*,*)
    write(*,'(A,F12.4,A)') ' Fin parameter (m):            ', m_val, ' 1/m'
    write(*,'(A,F12.4  )') ' mL parameter:                 ', mL
    write(*,'(A,F12.6,A)') ' Base Cross-Section Area (Ac): ', Ac * 1.0d6, ' mm2'
    write(*,'(A,F12.6,A)') ' Surface Area (Asurf):         ', As * 1.0d4, ' cm2'
    write(*,*)

    write(*,*) '================================================================'
    write(*,*) '                  PERFORMANCE RESULTS'
    write(*,*) '================================================================'
    write(*,*)
    write(*,'(A,F12.4,A)') ' Fin efficiency (eta_f):       ', eta_f * 100.0d0, ' %'
    write(*,'(A,F12.4,A)') ' Fin effectiveness (eps_f):    ', epsilon_f, ''
    write(*,'(A,F12.4,A)') ' Actual Heat Dissipated (Q):   ', Q_fin, ' W'
    write(*,'(A,F12.4,A)') ' Max Heat Dissipation (Qmax):  ', Qmax, ' W'
    write(*,'(A,F12.6  )') ' Biot number (Bi = h*Ac/P*k):  ', Biot_f
    write(*,*)

    write(*,*) '================================================================'
    write(*,*) '                 VALIDITY ASSESSMENT'
    write(*,*) '================================================================'
    write(*,*)
    if (Biot_f < 0.1d0) then
        write(*,'(A,F8.5,A)') ' [OK]    Bi = ', Biot_f, ' < 0.1  →  1-D fin assumption VALID.'
    else
        write(*,'(A,F8.5,A)') ' [WARN]  Bi = ', Biot_f, ' >= 0.1 →  1-D introduces error; use 2-D analysis.'
    end if
    write(*,*)
    if (eta_f > 0.9d0) then
        write(*,*) ' [EXCELLENT]  eta > 90% — highly efficient fin design.'
    else if (eta_f >= 0.7d0) then
        write(*,*) ' [GOOD]       eta 70–90% — acceptable performance.'
    else if (eta_f >= 0.5d0) then
        write(*,*) ' [MARGINAL]   eta 50–70% — consider increasing k or reducing L.'
    else
        write(*,*) ' [POOR]       eta < 50%  — fin is under-performing. Redesign recommended.'
    end if
    write(*,*)
    if (epsilon_f < 2.0d0) then
        write(*,*) ' [WARNING]    Effectiveness < 2 — fin usage may not be justified.'
    else if (epsilon_f > 10.0d0) then
        write(*,*) ' [EXCELLENT]  Effectiveness > 10 — fin greatly enhances heat transfer.'
    else
        write(*,*) ' [GOOD]       Effectiveness in range 2–10 — fin is beneficial.'
    end if
    write(*,*)

    ! ---------------------------------------------------------
    ! 4. TEMPERATURE PROFILE ALONG FIN
    ! ---------------------------------------------------------
    write(*,*) '================================================================'
    write(*,*) '                 TEMPERATURE PROFILE ALONG FIN'
    write(*,*) '================================================================'
    write(*,*) '  z [mm]    |  T [deg-C]    |  theta/theta_b'
    write(*,*) ' ---------------------------------------------------'

    n_points = 20
    do i = 0, n_points
        z_val = L_m * (real(i, 8) / real(n_points, 8)) ! base to tip coordinate
        x_val = L_m - z_val                            ! tip coordinate (L to 0)

        if (shape_id == 1 .or. shape_id == 4) then
            ! Uniform Cross section
            select case(tip_cond)
            case(1) ! Adiabatic
                theta_x = theta_b * (cosh(m_val * (L_m - z_val)) / cosh(mL))
            case(2) ! Convective
                theta_x = theta_b * (cosh(m_val * (L_m - z_val)) + &
                          (h_tip / (m_val * k_cond)) * sinh(m_val * (L_m - z_val))) / &
                          (cosh(mL) + (h_tip / (m_val * k_cond)) * sinh(mL))
            case(3) ! Prescribed T
                theta_L = T_tip - T_ambient
                theta_x = theta_b * ((theta_L / theta_b) * sinh(m_val * z_val) + &
                          sinh(m_val * (L_m - z_val))) / sinh(mL)
            case(4) ! Infinite
                theta_x = theta_b * exp(-m_val * z_val)
            end select
        else if (shape_id == 2) then
            ! Conical Spine
            if (x_val / L_m < 1.0d-6) then
                theta_x = theta_b * mL / I1_base
            else
                arg_z = 2.0d0 * m_val * sqrt(x_val * L_m)
                call bessel_I1(arg_z, I1_val)
                theta_x = theta_b * (x_val / L_m)**(-0.5d0) * (I1_val / I1_base)
            end if
        else if (shape_id == 3) then
            ! Parabolic Spine
            arg_z = (4.0d0 / 3.0d0) * m_val * L_m * (x_val / L_m)**1.5d0
            call bessel_I0(arg_z, I0_val)
            theta_x = theta_b * (I0_val / I0_base)
        end if

        T_x = T_ambient + theta_x
        write(*,'(2X,F10.3,2X,A,2X,F10.4,4X,A,2X,F8.5)') &
            z_val * 1000.0d0, '|', T_x, '|', theta_x / theta_b
    end do
    write(*,*)

    ! ---------------------------------------------------------
    ! 5. OPTIMIZATION AND RECOMMENDATIONS
    ! ---------------------------------------------------------
    write(*,*) '================================================================'
    write(*,*) '                OPTIMIZATION RECOMMENDATIONS'
    write(*,*) '================================================================'
    write(*,*)
    if (shape_id == 1 .or. shape_id == 4) then
        opt_L = 1.4192d0 / m_val
        write(*,'(A,F10.2,A)') ' Optimal length for max heat/volume:    ', opt_L * 1000.0d0, ' mm'
        write(*,'(A,F10.2,A)') ' Suggested change in length:            ', (opt_L - L_m) * 1000.0d0, ' mm'
    else if (shape_id == 2) then
        opt_L = 2.02d0 / m_val
        write(*,'(A,F10.2,A)') ' Optimal length for conical spine:      ', opt_L * 1000.0d0, ' mm'
        write(*,'(A,F10.2,A)') ' Suggested change in length:            ', (opt_L - L_m) * 1000.0d0, ' mm'
    else if (shape_id == 3) then
        opt_L = 2.12d0 / m_val
        write(*,'(A,F10.2,A)') ' Optimal length for parabolic spine:    ', opt_L * 1000.0d0, ' mm'
        write(*,'(A,F10.2,A)') ' Suggested change in length:            ', (opt_L - L_m) * 1000.0d0, ' mm'
    end if
    write(*,*)

    ! ---------------------------------------------------------
    ! 6. MATERIAL COMPARISON
    ! ---------------------------------------------------------
    write(*,*) '================================================================'
    write(*,*) '                MATERIAL PERFORMANCE COMPARISON'
    write(*,*) '================================================================'
    write(*,*) '  (Same geometry, h, base and ambient temperatures)'
    write(*,*)
    call compare_mat('Aluminum (6061-T6)')
    call compare_mat('Copper (pure)     ')
    call compare_mat('Steel (low-carbon)')
    call compare_mat('Brass (Cu65/Zn35) ')
    call compare_mat('Cast Iron         ')
    write(*,*)
    write(*,*) '================================================================'
    write(*,*) '               CALCULATION COMPLETE'
    write(*,*) '================================================================'

contains

    subroutine bessel_I0(x, res)
        implicit none
        real(8), intent(in)  :: x
        real(8), intent(out) :: res
        real(8) :: tx, p
        if (x <= 3.75d0) then
            tx = (x / 3.75d0)**2
            res = 1.0d0 + tx*(3.5156229d0 + tx*(3.0899424d0 + tx*(1.2067492d0 &
                + tx*(0.2659732d0 + tx*(0.0360768d0 + tx*0.0045813d0)))))
        else
            tx = 3.75d0 / x
            p = 0.39894228d0 + tx*(0.01328592d0 + tx*(0.00225319d0 &
                + tx*(-0.00157565d0 + tx*(0.00916281d0 + tx*(-0.02057706d0 &
                + tx*(0.02635537d0 + tx*(-0.01647633d0 + tx*0.00392377d0)))))))
            res = (exp(x) / sqrt(x)) * p
        end if
    end subroutine bessel_I0

    subroutine bessel_I1(x, res)
        implicit none
        real(8), intent(in)  :: x
        real(8), intent(out) :: res
        real(8) :: tx, p
        if (x <= 3.75d0) then
            tx = (x / 3.75d0)**2
            res = x * (0.5d0 + tx*(0.87890594d0 + tx*(0.51498869d0 &
                + tx*(0.15084934d0 + tx*(0.02658733d0 + tx*(0.00301532d0 &
                + tx*0.00032411d0))))))
        else
            tx = 3.75d0 / x
            p = 0.39894228d0 + tx*(-0.03988024d0 + tx*(-0.00362018d0 &
                + tx*(0.00163801d0 + tx*(-0.01031555d0 + tx*(0.02282967d0 &
                + tx*(-0.02895312d0 + tx*(0.01787654d0 - tx*0.00420059d0)))))))
            res = (exp(x) / sqrt(x)) * p
        end if
    end subroutine bessel_I1

    subroutine bessel_I2(x, res)
        implicit none
        real(8), intent(in)  :: x
        real(8), intent(out) :: res
        real(8) :: I0_val, I1_val
        if (x < 0.2d0) then
            res = (x**2 / 8.0d0) * (1.0d0 + (x**2 / 12.0d0) + (x**4 / 384.0d0))
        else
            call bessel_I0(x, I0_val)
            call bessel_I1(x, I1_val)
            res = I0_val - (2.0d0 / x) * I1_val
        end if
    end subroutine bessel_I2

    subroutine compare_mat(name)
        implicit none
        character(len=*), intent(in) :: name
        real(8) :: k_m, m_m, mL_m, Q_m, eta_m
        real(8) :: I0_bm, I1_bm, I2_bm, th_Lm

        ! 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 m
        if (shape_id == 1) then
            m_m = sqrt(4.0d0 * h_conv / (k_m * d_m))
        else if (shape_id == 4) then
            m_m = sqrt((h_conv * P) / (k_m * Ac))
        else
            m_m = sqrt(4.0d0 * h_conv / (k_m * d_m))
        end if

        mL_m = m_m * L_m

        if (shape_id == 1 .or. shape_id == 4) then
            select case(tip_cond)
            case(1) ! Adiabatic
                eta_m = tanh(mL_m) / mL_m
                Q_m = eta_m * h_conv * As * theta_b
            case(2) ! Convective
                Q_m = sqrt(h_conv * P * k_m * Ac) * theta_b * &
                      (sinh(mL_m) + (h_tip / (m_m * k_m)) * cosh(mL_m)) / &
                      (cosh(mL_m) + (h_tip / (m_m * k_m)) * sinh(mL_m))
                eta_m = Q_m / (h_conv * As * theta_b + h_tip * Ac * theta_b)
            case(3) ! Prescribed T
                th_Lm = T_tip - T_ambient
                Q_m = sqrt(h_conv * P * k_m * Ac) * theta_b * (cosh(mL_m) - th_Lm / theta_b) / sinh(mL_m)
                eta_m = Q_m / (h_conv * As * theta_b)
            case(4) ! Infinite
                Q_m = sqrt(h_conv * P * k_m * Ac) * theta_b
                eta_m = 1.0d0 / mL_m
            end select
        else if (shape_id == 2) then
            call bessel_I1(2.0d0 * mL_m, I1_bm)
            call bessel_I2(2.0d0 * mL_m, I2_bm)
            if (I1_bm > 1d-30) then
                eta_m = (2.0d0 / mL_m) * (I2_bm / I1_bm)
            else
                eta_m = 1.0d0
            end if
            Q_m = eta_m * h_conv * As * theta_b
        else if (shape_id == 3) then
            call bessel_I0((4.0d0 / 3.0d0) * mL_m, I0_bm)
            call bessel_I1((4.0d0 / 3.0d0) * mL_m, I1_bm)
            if (I0_bm > 1d-30) then
                eta_m = (1.5d0 / mL_m) * (I1_bm / I0_bm)
            else
                eta_m = 1.0d0
            end if
            Q_m = eta_m * h_conv * As * theta_b
        end if

        if (eta_m > 1.0d0) eta_m = 1.0d0
        if (eta_m < 0.0d0) eta_m = 0.0d0

        write(*,'(3X,A20,A,F6.1,A,F7.2,A,F9.3,A)') name, ': k = ', k_m, &
            ' W/m.K → eta = ', eta_m * 100.0d0, ' %, Q = ', Q_m, ' W'
    end subroutine compare_mat

end program Pin_Fin_Spine_Calculator