🔬
Solver Purpose & Physical Scope
Lumped system analysis simplifies transient heat transfer by assuming that the solid is at a uniform temperature at any instant during the process. This is valid when internal conduction resistance is much smaller than external convection resistance:
📂 Discipline: Conduction
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 243 times
📄 Source File:
lumped_system_analysis.f90
📁 calcul/Conduction /
lumped_system_analysis.f90
program lumped_system_analysis
implicit none
! Variable declarations
integer :: geometry_type, n_points, i
real(8) :: k, rho, cp, h, T_initial, T_ambient, T_target
real(8) :: radius, length_dim, side_dim
real(8) :: volume, A_surface, Lc, Bi, tau, alpha
real(8) :: T_t, t, dt, t_max, time_target, Q_t, Q_max
real(8) :: theta_ratio, pi
character(len=20) :: geometry_name
pi = 3.14159265358979d0
! =============================================
! READ INPUT PARAMETERS
! =============================================
read *, geometry_type ! 1=Sphere, 2=Long Cylinder, 3=Plane Wall, 4=Cube
select case (geometry_type)
case (1)
geometry_name = "SPHERE"
read *, radius ! Radius [m]
case (2)
geometry_name = "LONG CYLINDER"
read *, radius ! Radius [m]
read *, length_dim ! Length [m]
case (3)
geometry_name = "PLANE WALL"
read *, length_dim ! Half-thickness L [m]
read *, side_dim ! Area of one face [m2]
case (4)
geometry_name = "CUBE"
read *, side_dim ! Side length [m]
case default
print *, 'ERROR: Invalid geometry type'
stop
end select
read *, k ! Thermal conductivity [W/m.K]
read *, rho ! Density [kg/m3]
read *, cp ! Specific heat capacity [J/kg.K]
read *, h ! Convection coefficient [W/m2.K]
read *, T_initial ! Initial temperature [deg-C]
read *, T_ambient ! Ambient temperature [deg-C]
read *, T_target ! Target temperature [deg-C] (0 = no target)
! =============================================
! CALCULATE GEOMETRY
! =============================================
select case (geometry_type)
case (1) ! Sphere
volume = (4.0d0 / 3.0d0) * pi * radius**3
A_surface = 4.0d0 * pi * radius**2
case (2) ! Long Cylinder
volume = pi * radius**2 * length_dim
A_surface = 2.0d0 * pi * radius * length_dim + 2.0d0 * pi * radius**2
case (3) ! Plane Wall (2L thickness, area A on each face)
volume = 2.0d0 * length_dim * side_dim
A_surface = 2.0d0 * side_dim
case (4) ! Cube
volume = side_dim**3
A_surface = 6.0d0 * side_dim**2
end select
! Characteristic length
Lc = volume / A_surface
! Thermal diffusivity
alpha = k / (rho * cp)
! Biot number
Bi = h * Lc / k
! Time constant
tau = rho * cp * volume / (h * A_surface)
! Maximum energy transfer
Q_max = rho * cp * volume * abs(T_initial - T_ambient)
! =============================================
! DISPLAY RESULTS
! =============================================
print *, '================================================='
print *, ' LUMPED SYSTEM ANALYSIS'
print *, ' (Transient Heat Conduction)'
print *, '================================================='
print *, ''
print *, ' Geometry: ', trim(geometry_name)
print *, ''
print *, '================================================='
print *, ' INPUT PARAMETERS'
print *, '================================================='
print *, ''
select case (geometry_type)
case (1)
print '(A,F12.6,A)', ' Radius: ', radius, ' m'
print '(A,F12.4,A)', ' ', radius*1000.0d0, ' mm'
case (2)
print '(A,F12.6,A)', ' Radius: ', radius, ' m'
print '(A,F12.4,A)', ' Length: ', length_dim, ' m'
case (3)
print '(A,F12.6,A)', ' Half-thickness (L): ', length_dim, ' m'
print '(A,F12.4,A)', ' ', length_dim*1000.0d0, ' mm'
print '(A,F12.4,A)', ' Face Area (A): ', side_dim, ' m2'
case (4)
print '(A,F12.6,A)', ' Side Length: ', side_dim, ' m'
print '(A,F12.4,A)', ' ', side_dim*1000.0d0, ' mm'
end select
print '(A,F12.4,A)', ' Thermal Conductivity (k): ', k, ' W/m.K'
print '(A,F12.2,A)', ' Density (rho): ', rho, ' kg/m3'
print '(A,F12.2,A)', ' Specific Heat (cp): ', cp, ' J/kg.K'
print '(A,F12.2,A)', ' Convection Coeff. (h): ', h, ' W/m2.K'
print '(A,F12.2,A)', ' Initial Temperature (Ti): ', T_initial, ' deg-C'
print '(A,F12.2,A)', ' Ambient Temperature (T_inf):', T_ambient, ' deg-C'
if (T_target /= 0.0d0) then
print '(A,F12.2,A)', ' Target Temperature: ', T_target, ' deg-C'
end if
print *, ''
print *, '================================================='
print *, ' GEOMETRIC CALCULATIONS'
print *, '================================================='
print *, ''
print '(A,ES14.6,A)', ' Volume (V): ', volume, ' m3'
print '(A,ES14.6,A)', ' Surface Area (As): ', A_surface, ' m2'
print '(A,ES14.6,A)', ' Characteristic Length (Lc):', Lc, ' m'
print '(A,F12.4,A)', ' ', Lc*1000.0d0, ' mm'
print '(A,ES14.6,A)', ' Thermal Diffusivity (a): ', alpha, ' m2/s'
print *, ''
! =============================================
! BIOT NUMBER ANALYSIS
! =============================================
print *, '================================================='
print *, ' BIOT NUMBER ANALYSIS'
print *, '================================================='
print *, ''
print *, ' Bi = h * Lc / k'
print '(A,ES14.6)', ' Bi = ', Bi
print *, ''
if (Bi < 0.1d0) then
print *, ' STATUS: Bi < 0.1'
print *, ' -----------------------------------------------'
print *, ' The lumped system analysis IS VALID.'
print *, ' Internal temperature gradients are negligible.'
print *, ' The object can be treated as spatially uniform.'
else if (Bi < 1.0d0) then
print *, ' WARNING: 0.1 < Bi < 1.0'
print *, ' -----------------------------------------------'
print *, ' Lumped analysis has LIMITED accuracy.'
print *, ' Consider using the one-term approximation'
print *, ' (Heisler charts) for better precision.'
print *, ' Error may be up to 5%.'
else
print *, ' ERROR: Bi > 1.0'
print *, ' -----------------------------------------------'
print *, ' Lumped system analysis is NOT VALID!'
print *, ' Internal temperature gradients are significant.'
print *, ' Use the Heisler chart method or numerical'
print *, ' methods for accurate results.'
end if
print *, ''
! =============================================
! TIME CONSTANT
! =============================================
print *, '================================================='
print *, ' TIME CONSTANT'
print *, '================================================='
print *, ''
print *, ' tau = rho * cp * V / (h * As)'
print '(A,F14.4,A)', ' tau = ', tau, ' s'
if (tau >= 3600.0d0) then
print '(A,F14.4,A)', ' ', tau/3600.0d0, ' hours'
else if (tau >= 60.0d0) then
print '(A,F14.4,A)', ' ', tau/60.0d0, ' min'
end if
print *, ''
print *, ' Physical meaning:'
print *, ' At t = tau, the object reaches 63.2% of the'
print *, ' total temperature change.'
print *, ' At t = 3*tau, it reaches 95.0%.'
print *, ' At t = 5*tau, it reaches 99.3%.'
print *, ''
! =============================================
! TARGET TEMPERATURE
! =============================================
if (T_target /= 0.0d0) then
theta_ratio = (T_target - T_ambient) / (T_initial - T_ambient)
if (theta_ratio > 0.0d0 .and. theta_ratio < 1.0d0) then
time_target = -tau * log(theta_ratio)
print *, '================================================='
print *, ' TIME TO REACH TARGET TEMPERATURE'
print *, '================================================='
print *, ''
print '(A,F12.2,A)', ' Target Temperature: ', T_target, ' deg-C'
print '(A,F14.4,A)', ' Time Required: ', time_target, ' s'
if (time_target >= 3600.0d0) then
print '(A,F14.4,A)', ' ', time_target/3600.0d0, ' hours'
else if (time_target >= 60.0d0) then
print '(A,F14.4,A)', ' ', time_target/60.0d0, ' min'
end if
! Energy transferred at target time
Q_t = Q_max * (1.0d0 - exp(-time_target / tau))
print *, ''
print '(A,F14.4,A)', ' Energy Transferred: ', Q_t, ' J'
if (Q_t >= 1000.0d0) then
print '(A,F14.4,A)', ' ', Q_t/1000.0d0, ' kJ'
end if
print *, ''
else if (theta_ratio <= 0.0d0) then
print *, '================================================='
print *, ' TARGET TEMPERATURE ANALYSIS'
print *, '================================================='
print *, ''
print *, ' WARNING: Target temperature is beyond the'
print *, ' ambient temperature. The object will'
print *, ' asymptotically approach T_ambient but'
print *, ' never reach/cross it.'
print *, ''
else
print *, '================================================='
print *, ' TARGET TEMPERATURE ANALYSIS'
print *, '================================================='
print *, ''
print *, ' NOTE: Target temperature is between Ti'
print *, ' and T_ambient or equals Ti. Already reached.'
print *, ''
end if
end if
! =============================================
! ENERGY ANALYSIS
! =============================================
print *, '================================================='
print *, ' ENERGY ANALYSIS'
print *, '================================================='
print *, ''
print '(A,F14.4,A)', ' Max Energy Transfer (Qmax):', Q_max, ' J'
if (Q_max >= 1000.0d0) then
print '(A,F14.4,A)', ' ', Q_max/1000.0d0, ' kJ'
end if
if (Q_max >= 1.0d6) then
print '(A,F14.4,A)', ' ', Q_max/1.0d6, ' MJ'
end if
print *, ''
print *, ' Energy at key times:'
print '(A,F10.2,A,F14.4,A)', ' At t = 1*tau (', tau, ' s): Q = ', Q_max * 0.6321d0, ' J'
print '(A,F10.2,A,F14.4,A)', ' At t = 2*tau (', 2.0d0*tau, ' s): Q = ', Q_max * 0.8647d0, ' J'
print '(A,F10.2,A,F14.4,A)', ' At t = 3*tau (', 3.0d0*tau, ' s): Q = ', Q_max * 0.9502d0, ' J'
print '(A,F10.2,A,F14.4,A)', ' At t = 5*tau (', 5.0d0*tau, ' s): Q = ', Q_max * 0.9933d0, ' J'
print *, ''
! =============================================
! TEMPERATURE DISTRIBUTION OVER TIME
! =============================================
print *, '================================================='
print *, ' TEMPERATURE vs TIME'
print *, '================================================='
print *, ''
print *, ' Time (s) | T (deg-C) | theta/theta0 | Q/Qmax'
print *, ' ---------------------------------------------------------'
n_points = 50
t_max = 5.0d0 * tau
dt = t_max / dble(n_points)
do i = 0, n_points
t = dble(i) * dt
theta_ratio = exp(-t / tau)
T_t = T_ambient + (T_initial - T_ambient) * theta_ratio
Q_t = Q_max * (1.0d0 - theta_ratio)
write(*, '(3X,F12.4,4X,A,2X,F10.3,6X,A,2X,F10.6,4X,A,2X,F8.5)') &
t, '|', T_t, '|', theta_ratio, '|', Q_t / Q_max
end do
print *, ''
! =============================================
! HEAT TRANSFER RATE
! =============================================
print *, '================================================='
print *, ' INSTANTANEOUS HEAT TRANSFER RATE'
print *, '================================================='
print *, ''
print *, ' dQ/dt = h * As * (T(t) - T_inf)'
print *, ''
! At t=0
print '(A,F14.4,A)', ' At t = 0: dQ/dt = ', &
h * A_surface * abs(T_initial - T_ambient), ' W'
! At t=tau
print '(A,F14.4,A)', ' At t = tau: dQ/dt = ', &
h * A_surface * abs(T_initial - T_ambient) * exp(-1.0d0), ' W'
! At t=3tau
print '(A,F14.4,A)', ' At t = 3tau: dQ/dt = ', &
h * A_surface * abs(T_initial - T_ambient) * exp(-3.0d0), ' W'
print *, ''
! =============================================
! DESIGN RECOMMENDATIONS
! =============================================
print *, '================================================='
print *, ' DESIGN RECOMMENDATIONS'
print *, '================================================='
print *, ''
if (Bi < 0.1d0) then
print *, ' VALIDITY: Analysis is valid (Bi < 0.1)'
print *, ''
end if
! Suggestions to speed up or slow down
if (T_initial > T_ambient) then
print *, ' To SPEED UP cooling:'
else
print *, ' To SPEED UP heating:'
end if
print *, ' - Increase h (forced convection, fluid change)'
print *, ' - Decrease object size (reduce V/As ratio)'
print *, ' - Use material with lower rho*cp product'
print *, ''
if (T_initial > T_ambient) then
print *, ' To SLOW DOWN cooling:'
else
print *, ' To SLOW DOWN heating:'
end if
print *, ' - Add insulation (decrease h)'
print *, ' - Use larger objects (increase V/As ratio)'
print *, ' - Use material with higher rho*cp product'
print *, ''
! Convection regime recommendations
print *, ' Current convection regime:'
if (h < 10.0d0) then
print *, ' - Natural convection in air (h < 10)'
print *, ' - Consider forced air (h = 25-250) for faster results'
else if (h < 250.0d0) then
print *, ' - Forced convection in air (10 < h < 250)'
print *, ' - Consider liquid cooling (h > 500) if needed'
else if (h < 1000.0d0) then
print *, ' - Moderate liquid convection (250 < h < 1000)'
else
print *, ' - Strong liquid convection or boiling (h > 1000)'
end if
print *, ''
print *, '================================================='
print *, ' CALCULATION COMPLETE'
print *, '================================================='
print *, ''
end program lumped_system_analysis
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 lumped_system_analysis.f90 -o lumped_system
2. Execution with input.txt redirection:
lumped_system < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 0.005 63.9 7832.0 434.0 12.0 500.0 25.0 50.0
Parameter Description:
Geometry (1=Sphere, 2=Cylinder, 3=Plane Wall, 4=Cube) Radius/Half-thickness [m] Thermal conductivity k [W/m-K] Density rho [kg/m3] Specific heat Cp [J/kg-K] Convection coefficient h [W/m2-K] Initial temperature Ti [°C] Ambient temperature Tinf [°C] Target temperature T_target [°C]