💻 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.
Lumped System Analysis (Transient)
Core Numerical Engine in Fortran 90 • 0 total downloads
lumped_system_analysis.f90
! =========================================================================
! Source File: 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
Solver Description
Analyzes transient heat cooling/heating for solids with negligible internal temperature gradients. Verifies Biot number ($Bi = h L_c / k < 0.1$, where characteristic length $L_c = V/A_s$). Calculates temperature at time $t$ using the exponential time constant decay $\tau = \rho c_p V / (h A_s)$.
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 lumped_system_analysis.f90 -o lumped_system
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
lumped_system < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Geometry Code (1=Sphere, 2=Cylinder, 3=Plate)
1
! Characteristic Dimension (Radius or Half-Thickness) [m]
0.03
! Material Density ($\rho$) [kg/m³]
8500.0
! Specific Heat ($c_p$) [J/kg-K]
380.0
! Conductivity ($k$) [W/m-K]
110.0
! Convection Coefficient ($h$) [W/m²-K]
250.0
! Initial Temperature ($T_i$) [°C]
150.0
! Ambient Temperature ($T_\infty$) [°C]
20.0
! Time ($t$) [s]
180.0
1
! Characteristic Dimension (Radius or Half-Thickness) [m]
0.03
! Material Density ($\rho$) [kg/m³]
8500.0
! Specific Heat ($c_p$) [J/kg-K]
380.0
! Conductivity ($k$) [W/m-K]
110.0
! Convection Coefficient ($h$) [W/m²-K]
250.0
! Initial Temperature ($T_i$) [°C]
150.0
! Ambient Temperature ($T_\infty$) [°C]
20.0
! Time ($t$) [s]
180.0