๐ป Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Annular (Circular) Fin Efficiency
Core Numerical Engine in Fortran 90 โข 41 total downloads
annular_fin_efficiency.f90
! =========================================================================
! Source File: annular_fin_efficiency.f90
! =========================================================================
program Annular_Fin_Efficiency
implicit none
! ----------------------------------------------------------------
! ANNULAR (CIRCULAR) FIN EFFICIENCY CALCULATOR
!
! Correct formula from Incropera et al., "Fundamentals of Heat
! and Mass Transfer", 7th ed., Table 3.5, Case 4:
!
! 2 * r1 K1(m*r1) * I1(m*r2c) - I1(m*r1) * K1(m*r2c)
! eta_f = โโโโโโโโโโโ * โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
! m*(r2cยฒ-r1ยฒ) I0(m*r1) * K1(m*r2c) + K0(m*r1) * I1(m*r2c)
!
! where r2c = r2 + t/2 (corrected radius for convective tip)
! m = sqrt(2h / k*t)
!
! Modified Bessel functions computed via Abramowitz & Stegun ยง9.8
! ----------------------------------------------------------------
real(8) :: r1, t, r2, r2c, k, h, T0, Tinf
real(8) :: L, m, mR1, mR2c
real(8) :: I0_r1, I1_r1, K0_r1, K1_r1
real(8) :: I0_r2c, I1_r2c, K0_r2c, K1_r2c
real(8) :: num, den, eta
real(8) :: Asurf, Qmax, Q_fin, theta_base, epsilon_fin, Biot_fin
real(8) :: C1, C2, r_i, mr_i, I0_ri, K0_ri, theta_r, T_r
real(8), parameter :: PI = 3.14159265358979323846d0
integer :: i, n_points
write(*,*) '================================================================'
write(*,*) ' ANNULAR (CIRCULAR) FIN EFFICIENCY CALCULATOR'
write(*,*) '================================================================'
write(*,*)
write(*,*) 'Reading input parameters...'
! Inputs
read(*,*) r1 ! Inner radius [m]
read(*,*) t ! Fin thickness [m]
read(*,*) r2c ! Corrected tip radius r2c = r2 + t/2 [m]
read(*,*) k ! Thermal conductivity [W/m.K]
read(*,*) h ! Convection coefficient [W/m2.K]
read(*,*) T0 ! Base temperature [deg-C]
read(*,*) Tinf ! Ambient temperature [deg-C]
r2 = r2c - t / 2.0d0
L = r2 - r1
write(*,*)
write(*,*) '================================================================'
write(*,*) ' INPUT PARAMETERS'
write(*,*) '================================================================'
write(*,*)
write(*,'(A,F10.4,A)') ' Inner radius (r1): ', r1*1000.0d0, ' mm'
write(*,'(A,F10.4,A)') ' Outer radius (r2): ', r2*1000.0d0, ' mm'
write(*,'(A,F10.4,A)') ' Corrected tip radius (r2c): ', r2c*1000.0d0, ' mm'
write(*,'(A,F10.4,A)') ' Fin thickness (t): ', t*1000.0d0, ' mm'
write(*,'(A,F10.4,A)') ' Fin radial length (L=r2-r1): ', L*1000.0d0, ' mm'
write(*,'(A,F10.3,A)') ' Thermal Conductivity (k): ', k, ' W/m.K'
write(*,'(A,F10.3,A)') ' Convection Coefficient (h): ', h, ' W/m2.K'
write(*,'(A,F10.2,A)') ' Base Temperature (T0): ', T0, ' deg-C'
write(*,'(A,F10.2,A)') ' Ambient Temperature (Tinf): ', Tinf, ' deg-C'
write(*,*)
! Fin parameter
m = sqrt(2.0d0 * h / (k * t))
mR1 = m * r1
mR2c = m * r2c
write(*,*) '================================================================'
write(*,*) ' FIN PARAMETERS'
write(*,*) '================================================================'
write(*,*)
write(*,'(A,F12.4,A)') ' Fin parameter (m): ', m, ' 1/m'
write(*,'(A,F12.6 )') ' m * r1: ', mR1
write(*,'(A,F12.6 )') ' m * r2c (corrected): ', mR2c
write(*,'(A,F12.6 )') ' Radius ratio (r2/r1): ', r2 / r1
write(*,'(A,F12.6 )') ' Corrected ratio (r2c/r1): ', r2c / r1
write(*,*)
! Compute Bessel functions
call bessel_I0(mR1, I0_r1)
call bessel_I1(mR1, I1_r1)
call bessel_K0(mR1, K0_r1)
call bessel_K1(mR1, K1_r1)
call bessel_I0(mR2c, I0_r2c)
call bessel_I1(mR2c, I1_r2c)
call bessel_K0(mR2c, K0_r2c)
call bessel_K1(mR2c, K1_r2c)
write(*,*) '================================================================'
write(*,*) ' BESSEL FUNCTION VALUES'
write(*,*) '================================================================'
write(*,*)
write(*,'(A,F14.8)') ' I0(m*r1) = ', I0_r1
write(*,'(A,F14.8)') ' I1(m*r1) = ', I1_r1
write(*,'(A,F14.8)') ' K0(m*r1) = ', K0_r1
write(*,'(A,F14.8)') ' K1(m*r1) = ', K1_r1
write(*,*)
write(*,'(A,F14.8)') ' I0(m*r2c) = ', I0_r2c
write(*,'(A,F14.8)') ' I1(m*r2c) = ', I1_r2c
write(*,'(A,F14.8)') ' K0(m*r2c) = ', K0_r2c
write(*,'(A,F14.8)') ' K1(m*r2c) = ', K1_r2c
write(*,*)
! ----------------------------------------------------------------
! Fin efficiency (Incropera Table 3.5, Case 4):
!
! 2 r1 K1(mr1)*I1(mr2c) - I1(mr1)*K1(mr2c)
! eta_f = โโโโโโโโโ ร โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
! m(r2cยฒ-r1ยฒ) I0(mr1)*K1(mr2c) + K0(mr1)*I1(mr2c)
! ----------------------------------------------------------------
num = K1_r1 * I1_r2c - I1_r1 * K1_r2c
den = I0_r1 * K1_r2c + K0_r1 * I1_r2c
if (abs(den) < 1.0d-15) then
write(*,*) 'ERROR: Denominator near zero โ check geometry inputs.'
stop
end if
eta = (2.0d0 * r1 / (m * (r2c**2 - r1**2))) * (num / den)
if (eta > 1.0d0) eta = 1.0d0
if (eta < 0.0d0) eta = 0.0d0
! Performance
theta_base = T0 - Tinf
Asurf = 2.0d0 * PI * (r2c**2 - r1**2)
Qmax = h * Asurf * theta_base
Q_fin = eta * Qmax
epsilon_fin = Q_fin / max(h * 2.0d0 * PI * r1 * t * theta_base, 1.0d-30)
Biot_fin = h * t / (2.0d0 * k)
write(*,*) '================================================================'
write(*,*) ' PERFORMANCE RESULTS'
write(*,*) '================================================================'
write(*,*)
write(*,'(A,F12.4,A)') ' Fin efficiency (eta): ', eta * 100.0d0, ' %'
write(*,'(A,F12.6,A)') ' Fin corrected area (A_fin): ', Asurf * 1.0d4, ' cm2'
write(*,'(A,F12.4,A)') ' Max heat transfer (Q_max): ', Qmax, ' W'
write(*,'(A,F12.4,A)') ' Actual heat transfer (Q_fin): ', Q_fin, ' W'
write(*,'(A,F12.4 )') ' Fin effectiveness (epsilon): ', epsilon_fin
write(*,'(A,F12.6 )') ' Biot number (Bi = h*t/2k): ', Biot_fin
write(*,*)
! Validity
write(*,*) '================================================================'
write(*,*) ' VALIDITY ASSESSMENT'
write(*,*) '================================================================'
write(*,*)
if (Biot_fin < 0.1d0) then
write(*,'(A,F8.5,A)') ' [OK] Bi = ', Biot_fin, ' < 0.1 โ 1-D fin assumption VALID.'
else
write(*,'(A,F8.5,A)') ' [WARN] Bi = ', Biot_fin, ' >= 0.1 โ 1-D introduces error; use 2-D analysis.'
end if
write(*,*)
if (eta > 0.9d0) then
write(*,*) ' [EXCELLENT] eta > 90% โ highly efficient fin design.'
else if (eta >= 0.7d0) then
write(*,*) ' [GOOD] eta 70โ90% โ acceptable performance.'
else if (eta >= 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_fin < 2.0d0) then
write(*,*) ' [WARNING] Effectiveness < 2 โ fin usage may not be justified.'
else if (epsilon_fin > 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(*,*)
! Temperature profile
! theta(r)/theta_b = C1*I0(mr) + C2*K0(mr)
! BC1: theta(r1) = theta_b => C1*I0(mr1) + C2*K0(mr1) = 1
! BC2: dT/dr|r2c = 0 => C1*I1(mr2c) - C2*K1(mr2c) = 0
! => C2 = C1 * I1(mr2c) / K1(mr2c)
! => C1 = K1(mr2c) / [I0(mr1)*K1(mr2c) + I1(mr2c)*K0(mr1)]
C1 = K1_r2c / (I0_r1 * K1_r2c + I1_r2c * K0_r1)
C2 = C1 * I1_r2c / K1_r2c
write(*,*) '================================================================'
write(*,*) ' RADIAL TEMPERATURE DISTRIBUTION'
write(*,*) '================================================================'
write(*,*)
write(*,*) ' r [mm] | T [deg-C] | theta/theta_b'
write(*,*) ' ---------------------------------------------------'
n_points = 20
do i = 0, n_points
r_i = r1 + (real(i, 8) / real(n_points, 8)) * (r2 - r1)
mr_i = m * r_i
call bessel_I0(mr_i, I0_ri)
call bessel_K0(mr_i, K0_ri)
theta_r = C1 * I0_ri + C2 * K0_ri
T_r = Tinf + theta_r * theta_base
write(*,'(2X,F10.3,2X,A,2X,F10.4,4X,A,2X,F8.5)') &
r_i*1000.0d0, '|', T_r, '|', theta_r
end do
write(*,*)
! Material comparison
write(*,*) '================================================================'
write(*,*) ' MATERIAL COMPARISON'
write(*,*) '================================================================'
write(*,*)
write(*,*) ' (Same geometry, h, T0, Tinf โ material conductivity varies)'
write(*,*)
call material_eta('Aluminum (6061-T6) ', 200.0d0, m, k, r1, r2c)
call material_eta('Copper (pure) ', 385.0d0, m, k, r1, r2c)
call material_eta('Steel (low-carbon) ', 50.0d0, m, k, r1, r2c)
call material_eta('Brass (Cu65/Zn35) ', 110.0d0, m, k, r1, r2c)
call material_eta('Cast Iron ', 52.0d0, m, k, r1, r2c)
write(*,*)
write(*,*) '================================================================'
write(*,*) ' CALCULATION COMPLETE'
write(*,*) '================================================================'
contains
! Abramowitz & Stegun ยง9.8.1โ9.8.8 polynomial approximations
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_K0(x, res)
implicit none
real(8), intent(in) :: x
real(8), intent(out) :: res
real(8) :: tx, p, I0_val
if (x <= 2.0d0) then
tx = (x / 2.0d0)**2
call bessel_I0(x, I0_val)
p = -0.57721566d0 + tx*(0.42278420d0 + tx*(0.23069756d0 &
+ tx*(0.03488590d0 + tx*(0.00262698d0 + tx*(0.00010750d0 &
+ tx*0.0000074d0)))))
res = -log(x / 2.0d0) * I0_val + p
else
tx = 2.0d0 / x
p = 1.25331414d0 + tx*(-0.07832358d0 + tx*(0.02189568d0 &
+ tx*(-0.01062446d0 + tx*(0.00587872d0 + tx*(-0.00251540d0 &
+ tx*0.00053208d0)))))
res = (exp(-x) / sqrt(x)) * p
end if
end subroutine bessel_K0
subroutine bessel_K1(x, res)
implicit none
real(8), intent(in) :: x
real(8), intent(out) :: res
real(8) :: tx, p, I1_val
if (x <= 2.0d0) then
tx = (x / 2.0d0)**2
call bessel_I1(x, I1_val)
p = 1.0d0 + tx*(0.15443144d0 + tx*(-0.67278579d0 &
+ tx*(-0.18156897d0 + tx*(-0.01919402d0 + tx*(-0.00110404d0 &
- tx*0.00004686d0)))))
res = log(x / 2.0d0) * I1_val + p / x
else
tx = 2.0d0 / x
p = 1.25331414d0 + tx*(0.23498619d0 + tx*(-0.03655620d0 &
+ tx*(0.01504268d0 + tx*(-0.00780353d0 + tx*(0.00325614d0 &
- tx*0.00068245d0)))))
res = (exp(-x) / sqrt(x)) * p
end if
end subroutine bessel_K1
subroutine material_eta(name, k_mat, m_ref, k_ref, r1_in, r2c_in)
implicit none
character(len=*), intent(in) :: name
real(8), intent(in) :: k_mat, m_ref, k_ref, r1_in, r2c_in
real(8) :: m_mat, mR1m, mR2cm
real(8) :: I0r1, I1r1, K0r1, K1r1
real(8) :: I0r2c, I1r2c, K0r2c, K1r2c
real(8) :: num_m, den_m, eta_m
! Scale m: m = sqrt(2h/kt), so m_mat = m_ref * sqrt(k_ref / k_mat)
m_mat = m_ref * sqrt(k_ref / k_mat)
mR1m = m_mat * r1_in
mR2cm = m_mat * r2c_in
call bessel_I0(mR1m, I0r1)
call bessel_I1(mR1m, I1r1)
call bessel_K0(mR1m, K0r1)
call bessel_K1(mR1m, K1r1)
call bessel_I0(mR2cm, I0r2c)
call bessel_I1(mR2cm, I1r2c)
call bessel_K0(mR2cm, K0r2c)
call bessel_K1(mR2cm, K1r2c)
num_m = K1r1 * I1r2c - I1r1 * K1r2c
den_m = I0r1 * K1r2c + K0r1 * I1r2c
if (abs(den_m) < 1.0d-15) then
write(*,'(3X,A20,A)') name, ': near-singular โ skipped'
return
end if
eta_m = (2.0d0 * r1_in / (m_mat * (r2c_in**2 - r1_in**2))) * (num_m / den_m)
if (eta_m > 1.0d0) eta_m = 1.0d0
if (eta_m < 0.0d0) eta_m = 0.0d0
write(*,'(3X,A20,A,F7.2,A,F7.3,A)') name, ': k = ', k_mat, &
' W/m.K โ eta = ', eta_m * 100.0d0, ' %'
end subroutine material_eta
end program Annular_Fin_Efficiency
Solver Description
Calculate annular (circular) fin efficiency, effectiveness, and heat transfer rate using Bessel functions. Radial temperature profile and material comparison for heat sink design.
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
๐ ๏ธ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -O3 annular_fin_efficiency.f90 -o annular_fin
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
annular_fin < input.txt
๐ฅ Downloads & Local Files
Preview of the required input file (input.txt):
! Inner radius r1 [m]
0.025
! Fin thickness t [m]
0.002
! Corrected outer radius r2c [m]
0.065
! Thermal conductivity k [W/m-K]
200.0
! Convection coefficient h [W/m2-K]
25.0
! Base temperature T0 [รยฐC]
100.0
! Ambient temperature Tinf [รยฐC]
20.0
0.025
! Fin thickness t [m]
0.002
! Corrected outer radius r2c [m]
0.065
! Thermal conductivity k [W/m-K]
200.0
! Convection coefficient h [W/m2-K]
25.0
! Base temperature T0 [รยฐC]
100.0
! Ambient temperature Tinf [รยฐC]
20.0