🔬
Solver Purpose & Physical Scope
Steady-state conduction in solids with uniform volumetric heat generation ($\dot{q}$) is governed by the Poisson Equation. For 1D geometries, the temperature profile and center temperature $T_{max}$ are:
📂 Discipline: Conduction
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 499 times
📄 Source File:
heat_generation.f90
📁 calcul/Conduction /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 heat_generation.f90 -o heat_gen
2. Execution with input.txt redirection:
heat_gen < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 0.01 15.0 20000000.0 4000.0 25.0
Parameter Description:
Geometry (1=Plane Wall, 2=Cylinder, 3=Sphere) Characteristic dimension (L, ro) [m] Thermal conductivity k [W/m-K] Volumetric heat generation rate q_dot [W/m3] Convection coefficient h [W/m2-K] Ambient fluid temperature Tinf [°C]