💻 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.
Conduction with Internal Heat Generation
Core Numerical Engine in Fortran 90 • 0 total downloads
heat_generation.f90
! =========================================================================
! Source File: heat_generation.f90
! =========================================================================
program heat_generation
implicit none
integer :: geom_type, n_points, i
real(8) :: dim_char, k_mat, q_gen, h_conv, temp_inf
real(8) :: cyl_length
real(8) :: volume, area_surf, temp_surface, temp_max, delta_temp
real(8) :: q_total, x_pos, temp_x, ratio_x
real(8) :: pi
character(len=20) :: geom_name
pi = 3.14159265358979d0
! Read inputs
read *, geom_type ! 1=Plane Wall, 2=Long Cylinder, 3=Sphere
read *, dim_char ! Half-thickness L (wall) or radius r0 [m]
select case (geom_type)
case (1)
geom_name = "PLANE WALL"
case (2)
geom_name = "LONG CYLINDER"
read *, cyl_length ! Cylinder length [m]
case (3)
geom_name = "SPHERE"
case default
print *, 'ERROR: Invalid geometry type'
stop
end select
read *, k_mat ! Thermal conductivity [W/m.K]
read *, q_gen ! Volumetric heat generation [W/m3]
read *, h_conv ! Convection coefficient [W/m2.K]
read *, temp_inf ! Ambient temperature [deg-C]
! =============================================
! GEOMETRY CALCULATIONS
! =============================================
select case (geom_type)
case (1) ! Plane Wall (per unit area, symmetric)
volume = 2.0d0 * dim_char * 1.0d0 ! per m2 of wall face
area_surf = 2.0d0 * 1.0d0 ! both faces per m2
case (2) ! Long Cylinder
volume = pi * dim_char**2 * cyl_length
area_surf = 2.0d0 * pi * dim_char * cyl_length
case (3) ! Sphere
volume = (4.0d0 / 3.0d0) * pi * dim_char**3
area_surf = 4.0d0 * pi * dim_char**2
end select
! =============================================
! TEMPERATURE CALCULATIONS
! =============================================
select case (geom_type)
case (1) ! Plane Wall
temp_surface = temp_inf + q_gen * dim_char / h_conv
temp_max = temp_surface + q_gen * dim_char**2 / (2.0d0 * k_mat)
case (2) ! Long Cylinder
temp_surface = temp_inf + q_gen * dim_char / (2.0d0 * h_conv)
temp_max = temp_surface + q_gen * dim_char**2 / (4.0d0 * k_mat)
case (3) ! Sphere
temp_surface = temp_inf + q_gen * dim_char / (3.0d0 * h_conv)
temp_max = temp_surface + q_gen * dim_char**2 / (6.0d0 * k_mat)
end select
delta_temp = temp_max - temp_surface
q_total = q_gen * volume
! =============================================
! DISPLAY RESULTS
! =============================================
print *, '================================================='
print *, ' CONDUCTION WITH INTERNAL HEAT GENERATION'
print *, '================================================='
print *, ''
print *, ' Geometry: ', trim(geom_name)
print *, ''
print *, '================================================='
print *, ' INPUT PARAMETERS'
print *, '================================================='
print *, ''
select case (geom_type)
case (1)
print '(A,F12.6,A)', ' Half-thickness (L): ', dim_char, ' m'
print '(A,F12.4,A)', ' ', dim_char*1000.0d0, ' mm'
case (2)
print '(A,F12.6,A)', ' Radius (r0): ', dim_char, ' m'
print '(A,F12.4,A)', ' ', dim_char*1000.0d0, ' mm'
print '(A,F12.4,A)', ' Cylinder Length: ', cyl_length, ' m'
case (3)
print '(A,F12.6,A)', ' Radius (r0): ', dim_char, ' m'
print '(A,F12.4,A)', ' ', dim_char*1000.0d0, ' mm'
end select
print '(A,F12.4,A)', ' Thermal Conductivity (k): ', k_mat, ' W/m.K'
print '(A,ES14.4,A)', ' Heat Generation (q_gen): ', q_gen, ' W/m3'
print '(A,F12.2,A)', ' Convection Coeff. (h): ', h_conv, ' W/m2.K'
print '(A,F12.2,A)', ' Ambient Temp. (T_inf): ', temp_inf, ' deg-C'
print *, ''
print *, '================================================='
print *, ' GEOMETRIC PROPERTIES'
print *, '================================================='
print *, ''
print '(A,ES14.6,A)', ' Volume: ', volume, ' m3'
print '(A,ES14.6,A)', ' Surface Area: ', area_surf, ' m2'
print *, ''
! =============================================
! TEMPERATURE RESULTS
! =============================================
print *, '================================================='
print *, ' TEMPERATURE DISTRIBUTION'
print *, '================================================='
print *, ''
select case (geom_type)
case (1)
print *, ' Formula (Plane Wall, symmetric):'
print *, ' T(x) = Ts + q_gen/(2k) * (L^2 - x^2)'
print *, ' Ts = T_inf + q_gen*L/h'
case (2)
print *, ' Formula (Long Cylinder):'
print *, ' T(r) = Ts + q_gen/(4k) * (r0^2 - r^2)'
print *, ' Ts = T_inf + q_gen*r0/(2h)'
case (3)
print *, ' Formula (Sphere):'
print *, ' T(r) = Ts + q_gen/(6k) * (r0^2 - r^2)'
print *, ' Ts = T_inf + q_gen*r0/(3h)'
end select
print *, ''
print '(A,F12.2,A)', ' Surface Temperature (Ts): ', temp_surface, ' deg-C'
print '(A,F12.2,A)', ' Maximum Temperature (Tmax):', temp_max, ' deg-C'
print '(A,F12.2,A)', ' Temperature Drop (dT): ', delta_temp, ' deg-C'
print '(A,F12.2,A)', ' (Tmax - Ts)'
print *, ''
! =============================================
! HEAT TRANSFER
! =============================================
print *, '================================================='
print *, ' HEAT TRANSFER ANALYSIS'
print *, '================================================='
print *, ''
print '(A,ES14.6,A)', ' Total Heat Generated: ', q_total, ' W'
if (q_total >= 1000.0d0) then
print '(A,F14.4,A)', ' ', q_total / 1000.0d0, ' kW'
end if
print '(A,ES14.6,A)', ' Heat Flux at Surface: ', q_total / area_surf, ' W/m2'
print *, ''
! Energy balance verification
print *, ' Energy Balance Verification:'
print '(A,ES14.6,A)', ' Q_generated = q_gen * V = ', q_total, ' W'
print '(A,ES14.6,A)', ' Q_convection = h*As*dTs = ', &
h_conv * area_surf * (temp_surface - temp_inf), ' W'
print *, ''
! =============================================
! TEMPERATURE PROFILE
! =============================================
print *, '================================================='
print *, ' TEMPERATURE PROFILE'
print *, '================================================='
print *, ''
select case (geom_type)
case (1)
print *, ' x/L | x (mm) | T (deg-C) | (T-Ts)/dT'
case default
print *, ' r/r0 | r (mm) | T (deg-C) | (T-Ts)/dT'
end select
print *, ' --------------------------------------------------------'
n_points = 25
do i = 0, n_points
ratio_x = dble(i) / dble(n_points)
x_pos = ratio_x * dim_char
select case (geom_type)
case (1)
temp_x = temp_surface + q_gen / (2.0d0 * k_mat) * &
(dim_char**2 - x_pos**2)
case (2)
temp_x = temp_surface + q_gen / (4.0d0 * k_mat) * &
(dim_char**2 - x_pos**2)
case (3)
temp_x = temp_surface + q_gen / (6.0d0 * k_mat) * &
(dim_char**2 - x_pos**2)
end select
write(*, '(3X,F8.4,5X,A,2X,F10.4,5X,A,2X,F10.2,5X,A,2X,F8.5)') &
ratio_x, '|', x_pos*1000.0d0, '|', temp_x, '|', &
(temp_x - temp_surface) / max(delta_temp, 1.0d-10)
end do
print *, ''
! =============================================
! THERMAL RESISTANCE ANALYSIS
! =============================================
print *, '================================================='
print *, ' THERMAL RESISTANCE BREAKDOWN'
print *, '================================================='
print *, ''
select case (geom_type)
case (1)
print '(A,ES14.6,A)', ' R_conduction (L/k): ', &
dim_char / k_mat, ' m2.K/W'
print '(A,ES14.6,A)', ' R_convection (1/h): ', &
1.0d0 / h_conv, ' m2.K/W'
print '(A,F12.4)', ' Ratio R_cond/R_conv: ', &
(dim_char / k_mat) / (1.0d0 / h_conv)
case (2)
print '(A,ES14.6,A)', ' R_convection: ', &
1.0d0 / (h_conv * area_surf), ' K/W'
case (3)
print '(A,ES14.6,A)', ' R_convection: ', &
1.0d0 / (h_conv * area_surf), ' K/W'
end select
print *, ''
! =============================================
! SAFETY ANALYSIS
! =============================================
print *, '================================================='
print *, ' SAFETY ANALYSIS'
print *, '================================================='
print *, ''
print '(A,F12.2,A)', ' Maximum Temperature: ', temp_max, ' deg-C'
print *, ''
if (temp_max > 1000.0d0) then
print *, ' CRITICAL: Tmax > 1000 deg-C'
print *, ' -------------------------------------------'
print *, ' Temperature exceeds most material limits!'
print *, ' Risk of melting or structural failure.'
else if (temp_max > 500.0d0) then
print *, ' WARNING: Tmax > 500 deg-C'
print *, ' -------------------------------------------'
print *, ' High temperature may cause material'
print *, ' degradation. Verify material compatibility.'
else if (temp_max > 200.0d0) then
print *, ' CAUTION: Tmax > 200 deg-C'
print *, ' -------------------------------------------'
print *, ' Moderate temperature. Check thermal limits'
print *, ' for polymers and organic materials.'
else
print *, ' STATUS: Tmax < 200 deg-C'
print *, ' -------------------------------------------'
print *, ' Temperature is within safe limits for'
print *, ' most engineering materials.'
end if
print *, ''
! =============================================
! PARAMETRIC STUDY
! =============================================
print *, '================================================='
print *, ' PARAMETRIC SENSITIVITY'
print *, '================================================='
print *, ''
print *, ' Effect of doubling key parameters on Tmax:'
print *, ''
! Double q_gen
select case (geom_type)
case (1)
print '(A,F12.2,A)', ' 2x q_gen -> Tmax = ', &
temp_inf + 2.0d0*q_gen*dim_char/h_conv + &
2.0d0*q_gen*dim_char**2/(2.0d0*k_mat), ' deg-C'
case (2)
print '(A,F12.2,A)', ' 2x q_gen -> Tmax = ', &
temp_inf + 2.0d0*q_gen*dim_char/(2.0d0*h_conv) + &
2.0d0*q_gen*dim_char**2/(4.0d0*k_mat), ' deg-C'
case (3)
print '(A,F12.2,A)', ' 2x q_gen -> Tmax = ', &
temp_inf + 2.0d0*q_gen*dim_char/(3.0d0*h_conv) + &
2.0d0*q_gen*dim_char**2/(6.0d0*k_mat), ' deg-C'
end select
! Double h
select case (geom_type)
case (1)
print '(A,F12.2,A)', ' 2x h -> Tmax = ', &
temp_inf + q_gen*dim_char/(2.0d0*h_conv) + &
q_gen*dim_char**2/(2.0d0*k_mat), ' deg-C'
case (2)
print '(A,F12.2,A)', ' 2x h -> Tmax = ', &
temp_inf + q_gen*dim_char/(4.0d0*h_conv) + &
q_gen*dim_char**2/(4.0d0*k_mat), ' deg-C'
case (3)
print '(A,F12.2,A)', ' 2x h -> Tmax = ', &
temp_inf + q_gen*dim_char/(6.0d0*h_conv) + &
q_gen*dim_char**2/(6.0d0*k_mat), ' deg-C'
end select
! Double k
select case (geom_type)
case (1)
print '(A,F12.2,A)', ' 2x k -> Tmax = ', &
temp_inf + q_gen*dim_char/h_conv + &
q_gen*dim_char**2/(4.0d0*k_mat), ' deg-C'
case (2)
print '(A,F12.2,A)', ' 2x k -> Tmax = ', &
temp_inf + q_gen*dim_char/(2.0d0*h_conv) + &
q_gen*dim_char**2/(8.0d0*k_mat), ' deg-C'
case (3)
print '(A,F12.2,A)', ' 2x k -> Tmax = ', &
temp_inf + q_gen*dim_char/(3.0d0*h_conv) + &
q_gen*dim_char**2/(12.0d0*k_mat), ' deg-C'
end select
print *, ''
print '(A,F12.2,A)', ' Current Tmax = ', temp_max, ' deg-C'
print *, ''
! =============================================
! RECOMMENDATIONS
! =============================================
print *, '================================================='
print *, ' DESIGN RECOMMENDATIONS'
print *, '================================================='
print *, ''
print *, ' To REDUCE maximum temperature:'
print *, ' - Decrease object size (reduce L or r0)'
print *, ' - Increase thermal conductivity (k)'
print *, ' - Increase convection coefficient (h)'
print *, ' - Reduce heat generation rate (q_gen)'
print *, ''
if (delta_temp > (temp_surface - temp_inf)) then
print *, ' NOTE: Internal resistance dominates'
print *, ' -> Focus on increasing k or reducing size'
else
print *, ' NOTE: External resistance dominates'
print *, ' -> Focus on increasing h (better cooling)'
end if
print *, ''
print *, '================================================='
print *, ' CALCULATION COMPLETE'
print *, '================================================='
print *, ''
end program heat_generation
Solver Description
Solves steady-state 1D heat conduction with uniform internal energy generation. Computes maximum center temperature ($T_{max}$), surface temperature ($T_s$), and local temperature distribution based on geometry-specific Poisson equation solutions.
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_generation.f90 -o heat_gen
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
heat_gen < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Geometry Code (1=Plane wall, 2=Cylinder, 3=Sphere)
1
! Half-Thickness or Outer Radius ($L$) [m]
0.02
! Thermal Conductivity ($k$) [W/m-K]
1.5
! Volumetric Heat Generation Rate (q_dot) [W/m³]
200000.0
! Convection Coefficient ($h$) [W/m²-K]
500.0
! Ambient Temperature ($T_\infty$) [°C]
25.0
1
! Half-Thickness or Outer Radius ($L$) [m]
0.02
! Thermal Conductivity ($k$) [W/m-K]
1.5
! Volumetric Heat Generation Rate (q_dot) [W/m³]
200000.0
! Convection Coefficient ($h$) [W/m²-K]
500.0
! Ambient Temperature ($T_\infty$) [°C]
25.0