💻 Fortran Source Code Library
Run calculations locally on your own machine. View code structure, read technical explanations, and download compilation packages including sample input files.
Composite Cylinder/Sphere Conduction
Core Numerical Engine in Fortran 90 • 0 total downloads
heat_conduction_cylinders_spheres.f90
! =========================================================================
! Source File: heat_conduction_cylinders_spheres.f90
! =========================================================================
program steady_heat_conduction_cylinders_spheres
implicit none
! Variable declarations
integer :: n_layers, i, n_nodes, j, geometry_type
real(8) :: T_inner, T_outer, Q, R_total, T_interface
real(8), allocatable :: r_inner(:), r_outer(:), k(:), R(:), T(:)
real(8) :: r_val, dr_val, length_cyl, heat_flux, r_mid
real(8) :: surface_inner, surface_outer
character(len=50) :: filename
character(len=20) :: geometry_name
! Read geometry type: 1 = Cylinder, 2 = Sphere
read *, geometry_type
if (geometry_type == 1) then
geometry_name = "CYLINDER"
read *, length_cyl ! Length for cylinder
else
geometry_name = "SPHERE"
length_cyl = 1.0d0 ! Not used for sphere
end if
! Read number of layers
read *, n_layers
! Allocate arrays
allocate(r_inner(n_layers))
allocate(r_outer(n_layers))
allocate(k(n_layers))
allocate(R(n_layers))
allocate(T(0:n_layers))
! Read layer properties (inner radius, outer radius, conductivity)
do i = 1, n_layers
read *, r_inner(i)
read *, r_outer(i)
read *, k(i)
end do
! Read boundary temperatures
read *, T_inner
read *, T_outer
! Calculate thermal resistances
R_total = 0.0d0
do i = 1, n_layers
if (geometry_type == 1) then
! Cylindrical resistance: R = ln(r_o/r_i) / (2*pi*k*L)
R(i) = log(r_outer(i) / r_inner(i)) / (2.0d0 * 3.14159265358979d0 * k(i) * length_cyl)
else
! Spherical resistance: R = (1/r_i - 1/r_o) / (4*pi*k)
R(i) = (1.0d0/r_inner(i) - 1.0d0/r_outer(i)) / (4.0d0 * 3.14159265358979d0 * k(i))
end if
R_total = R_total + R(i)
end do
! Calculate heat transfer rate
Q = (T_inner - T_outer) / R_total
! Calculate interface temperatures
T(0) = T_inner
T(n_layers) = T_outer
do i = 1, n_layers - 1
T(i) = T_inner - (T_inner - T_outer) * sum(R(1:i)) / R_total
end do
! Calculate surface areas
if (geometry_type == 1) then
surface_inner = 2.0d0 * 3.14159265358979d0 * r_inner(1) * length_cyl
surface_outer = 2.0d0 * 3.14159265358979d0 * r_outer(n_layers) * length_cyl
else
surface_inner = 4.0d0 * 3.14159265358979d0 * r_inner(1)**2
surface_outer = 4.0d0 * 3.14159265358979d0 * r_outer(n_layers)**2
end if
! Display results
print *, '========================================='
print *, 'STEADY STATE THERMAL CONDUCTION'
print *, 'GEOMETRY: ', trim(geometry_name)
print *, '========================================='
print *, ''
print *, 'INPUT PARAMETERS:'
print *, '----------------------------------------'
if (geometry_type == 1) then
print '(A,F10.4,A)', ' Cylinder Length: ', length_cyl, ' m'
print '(A,F10.4,A)', ' Inner Radius: ', r_inner(1), ' m'
print '(A,F10.4,A)', ' Outer Radius: ', r_outer(n_layers), ' m'
else
print '(A,F10.4,A)', ' Inner Radius: ', r_inner(1), ' m'
print '(A,F10.4,A)', ' Outer Radius: ', r_outer(n_layers), ' m'
end if
print '(A,F10.2,A)', ' Inner Temperature: ', T_inner, ' °C'
print '(A,F10.2,A)', ' Outer Temperature: ', T_outer, ' °C'
print '(A,F10.2,A)', ' Temperature Difference: ', (T_inner - T_outer), ' °C'
print '(A,I3)', ' Number of Layers: ', n_layers
print *, ''
print *, '========================================='
print *, 'MAIN RESULTS'
print *, '========================================='
print *, ''
print '(A,F12.6,A)', ' Total Thermal Resistance: ', R_total, ' K/W'
print '(A,F12.2,A)', ' Total Heat Flux (Q): ', Q, ' W'
print '(A,F12.4,A)', ' Inner Surface Area: ', surface_inner, ' m²'
print '(A,F12.4,A)', ' Outer Surface Area: ', surface_outer, ' m²'
print '(A,F12.2,A)', ' Inner Surface Heat Flux: ', Q/surface_inner, ' W/m²'
print '(A,F12.2,A)', ' Outer Surface Heat Flux: ', Q/surface_outer, ' W/m²'
print *, ''
print *, '========================================='
print *, 'LAYER ANALYSIS'
print *, '========================================='
print *, ''
do i = 1, n_layers
r_mid = (r_inner(i) + r_outer(i)) / 2.0d0
print '(A,I2,A)', ' === Layer ', i, ' ==='
print '(A,F10.4,A)', ' Inner Radius (r_i): ', r_inner(i), ' m'
print '(A,F10.4,A)', ' ', r_inner(i)*1000, ' mm'
print '(A,F10.4,A)', ' Outer Radius (r_o): ', r_outer(i), ' m'
print '(A,F10.4,A)', ' ', r_outer(i)*1000, ' mm'
print '(A,F10.4,A)', ' Thickness (r_o - r_i): ', (r_outer(i)-r_inner(i)), ' m'
print '(A,F10.4,A)', ' ', (r_outer(i)-r_inner(i))*1000, ' mm'
print '(A,F10.4,A)', ' Conductivity (k): ', k(i), ' W/m·K'
print '(A,F10.6,A)', ' Resistance (R): ', R(i), ' K/W'
print '(A,F8.2,A)', ' % of Total Resistance: ', (R(i)/R_total)*100, ' %'
print '(A,F10.2,A)', ' Inner Temperature: ', T(i-1), ' °C'
print '(A,F10.2,A)', ' Outer Temperature: ', T(i), ' °C'
print '(A,F10.2,A)', ' Temperature Drop: ', (T(i-1) - T(i)), ' °C'
print *, ''
end do
print *, '========================================='
print *, 'INTERFACE TEMPERATURES'
print *, '========================================='
print *, ''
print '(A,F10.2,A)', ' Inner Surface (T0): ', T(0), ' °C'
do i = 1, n_layers - 1
print '(A,I2,A,I2,A,F10.2,A)', ' Interface layer ', i, '-', i+1, ': ', T(i), ' °C'
end do
print '(A,F10.2,A)', ' Outer Surface (Tn): ', T(n_layers), ' °C'
print *, ''
! Generate temperature profile data file
filename = 'temperature_profile.dat'
open(unit=10, file=filename, status='replace')
write(10, '(A)') '# Radius(m) Temperature(°C) Layer'
do i = 1, n_layers
n_nodes = 50
dr_val = (r_outer(i) - r_inner(i)) / real(n_nodes, 8)
do j = 0, n_nodes
r_val = r_inner(i) + j * dr_val
if (geometry_type == 1) then
! Cylindrical temperature distribution
if (abs(r_outer(i) - r_inner(i)) > 1.0d-10) then
T_interface = T(i-1) - (T(i-1) - T(i)) * &
log(r_val/r_inner(i)) / log(r_outer(i)/r_inner(i))
else
T_interface = T(i-1)
end if
else
! Spherical temperature distribution
if (abs(r_outer(i) - r_inner(i)) > 1.0d-10) then
T_interface = T(i-1) - (T(i-1) - T(i)) * &
(1.0d0/r_inner(i) - 1.0d0/r_val) / &
(1.0d0/r_inner(i) - 1.0d0/r_outer(i))
else
T_interface = T(i-1)
end if
end if
write(10, '(F12.6,2X,F12.4,2X,I3)') r_val, T_interface, i
end do
end do
close(10)
print *, '========================================='
print *, 'DATA FILE'
print *, '========================================='
print *, ''
print *, ' Temperature profile saved in:'
print *, ' ', trim(filename)
print *, ''
print *, ' Format: Radius(m) Temperature(°C) Layer'
print *, ''
! Additional performance metrics
print *, '========================================='
print *, 'PERFORMANCE METRICS'
print *, '========================================='
print *, ''
! Find the layer with maximum resistance (bottleneck)
i = maxloc(R, dim=1)
print '(A,I2)', ' Limiting Layer (max resistance): Layer ', i
print '(A,F10.6,A)', ' Resistance of this layer: ', R(i), ' K/W'
print *, ''
! Calculate overall heat transfer coefficient
print '(A,F10.4,A)', ' U-Coefficient (inner basis): ', 1.0d0/(R_total*surface_inner), ' W/m²·K'
print '(A,F10.4,A)', ' U-Coefficient (outer basis): ', 1.0d0/(R_total*surface_outer), ' W/m²·K'
print *, ''
! Energy loss calculations
print *, 'ENERGY LOSSES (estimates):'
print '(A,F12.2,A)', ' Per hour: ', Q * 3600.0d0 / 1000.0d0, ' kJ/h'
print '(A,F12.2,A)', ' Per day: ', Q * 86400.0d0 / 1000.0d0, ' kJ/day'
print '(A,F12.2,A)', ' ', Q * 24.0d0 / 1000.0d0, ' kWh/day'
print '(A,F12.2,A)', ' Per year: ', Q * 8760.0d0 / 1000.0d0, ' kWh/year'
print *, ''
! Volume calculations
if (geometry_type == 1) then
print *, 'VOLUMETRIC CHARACTERISTICS:'
print '(A,F12.6,A)', ' Inner Volume: ', &
3.14159265358979d0 * r_inner(1)**2 * length_cyl, ' m³'
print '(A,F12.6,A)', ' Outer Volume: ', &
3.14159265358979d0 * r_outer(n_layers)**2 * length_cyl, ' m³'
print '(A,F12.6,A)', ' Wall Volume: ', &
3.14159265358979d0 * (r_outer(n_layers)**2 - r_inner(1)**2) * length_cyl, ' m³'
else
print *, 'VOLUMETRIC CHARACTERISTICS:'
print '(A,F12.6,A)', ' Inner Volume: ', &
(4.0d0/3.0d0) * 3.14159265358979d0 * r_inner(1)**3, ' m³'
print '(A,F12.6,A)', ' Outer Volume: ', &
(4.0d0/3.0d0) * 3.14159265358979d0 * r_outer(n_layers)**3, ' m³'
print '(A,F12.6,A)', ' Wall Volume: ', &
(4.0d0/3.0d0) * 3.14159265358979d0 * &
(r_outer(n_layers)**3 - r_inner(1)**3), ' m³'
end if
print *, ''
print *, '========================================='
print *, 'END OF CALCULATION'
print *, '========================================='
! Cleanup
deallocate(r_inner, r_outer, k, R, T)
end program steady_heat_conduction_cylinders_spheres
Solver Description
Solves steady-state 1D heat transfer in radial composite shells (cylindrical pipes or spherical containers). Sums radial resistances logarithmically for cylinders ($R = \ln(r_{out}/r_{in}) / (2\pi k L)$) or algebraically for spheres ($R = (r_{out}-r_{in}) / (4\pi k r_{in} r_{out})$) to compute heat loss rate ($Q$) and temperature profiles.
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 -ffree-line-length-none heat_conduction_cylinders_spheres.f90 -o heat_conduction_cyl_sph
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
heat_conduction_cyl_sph < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Geometry Code (1=Cylinder, 2=Sphere)
1
! Length ($L$) [m] (ignored for Sphere)
1.5
! Inner Radius ($r_i$) [m]
0.05
! Outer Radius ($r_o$) [m]
0.10
! Inner Surface Temp ($T_1$) [°C]
120.0
! Outer Surface Temp ($T_2$) [°C]
40.0
! Number of layers (N)
1
! Layer 1 Outer Radius ($r_1$) [m]
0.10
! Layer 1 Conductivity ($k_1$) [W/m-K]
45.0
1
! Length ($L$) [m] (ignored for Sphere)
1.5
! Inner Radius ($r_i$) [m]
0.05
! Outer Radius ($r_o$) [m]
0.10
! Inner Surface Temp ($T_1$) [°C]
120.0
! Outer Surface Temp ($T_2$) [°C]
40.0
! Number of layers (N)
1
! Layer 1 Outer Radius ($r_1$) [m]
0.10
! Layer 1 Conductivity ($k_1$) [W/m-K]
45.0