💻 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.
1D Transient Conduction (Heisler Series)
Core Numerical Engine in Fortran 90 • 0 total downloads
transient_conduction_1d.f90
! =========================================================================
! Source File: transient_conduction_1d.f90
! =========================================================================
program transient_conduction_1d
implicit none
! Variable declarations
integer :: geom_type, n_points, i, n_terms
real(8) :: dim_char, k_mat, rho_mat, cp_mat, h_conv
real(8) :: temp_init, temp_inf
real(8) :: alpha_diff, biot, fourier, time_val
real(8) :: x_ratio
real(8) :: zeta1, coeff_c1, theta_star, theta_center
real(8) :: temp_result, q_ratio
real(8) :: pi, fo_val, x_pos, temp_x
real(8) :: zeta_n, delta_fo
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
read *, k_mat ! Thermal conductivity [W/m.K]
read *, rho_mat ! Density [kg/m3]
read *, cp_mat ! Specific heat [J/kg.K]
read *, h_conv ! Convection coefficient [W/m2.K]
read *, temp_init ! Initial temperature [deg-C]
read *, temp_inf ! Ambient temperature [deg-C]
read *, x_ratio ! Position ratio x/L or r/r0 (0=center, 1=surface)
read *, time_val ! Time [s]
select case (geom_type)
case (1)
geom_name = "PLANE WALL"
case (2)
geom_name = "LONG CYLINDER"
case (3)
geom_name = "SPHERE"
case default
print *, 'ERROR: Invalid geometry type'
stop
end select
! Thermal diffusivity
alpha_diff = k_mat / (rho_mat * cp_mat)
! Biot number
biot = h_conv * dim_char / k_mat
! Fourier number
if (time_val > 0.0d0) then
fourier = alpha_diff * time_val / (dim_char**2)
else
fourier = 0.0d0
end if
! =============================================
! DISPLAY HEADER
! =============================================
print *, '================================================='
print *, ' TRANSIENT CONDUCTION - ONE-TERM APPROXIMATION'
print *, ' (Heisler Chart Method)'
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'
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,F12.2,A)', ' Density (rho): ', rho_mat, ' kg/m3'
print '(A,F12.2,A)', ' Specific Heat (cp): ', cp_mat, ' J/kg.K'
print '(A,F12.2,A)', ' Convection Coeff. (h): ', h_conv, ' W/m2.K'
print '(A,F12.2,A)', ' Initial Temp. (Ti): ', temp_init, ' deg-C'
print '(A,F12.2,A)', ' Ambient Temp. (T_inf): ', temp_inf, ' deg-C'
print '(A,F12.4)', ' Position ratio (x/L or r/r0):', x_ratio
print '(A,F12.4,A)', ' Time: ', time_val, ' s'
print *, ''
! =============================================
! DIMENSIONLESS NUMBERS
! =============================================
print *, '================================================='
print *, ' DIMENSIONLESS PARAMETERS'
print *, '================================================='
print *, ''
print '(A,ES14.6,A)', ' Thermal Diffusivity (a): ', alpha_diff, ' m2/s'
print '(A,F14.6)', ' Biot Number (Bi): ', biot
print '(A,F14.6)', ' Fourier Number (Fo): ', fourier
print *, ''
! Validity check
if (fourier < 0.2d0 .and. time_val > 0.0d0) then
print *, ' WARNING: Fo < 0.2'
print *, ' -------------------------------------------'
print *, ' The one-term approximation may not be'
print *, ' accurate. Consider using the full series'
print *, ' solution or numerical methods.'
print *, ' Accuracy improves significantly for Fo > 0.2'
else if (time_val > 0.0d0) then
print *, ' STATUS: Fo >= 0.2'
print *, ' -------------------------------------------'
print *, ' One-term approximation IS VALID.'
print *, ' Error is typically less than 2%.'
end if
print *, ''
if (biot < 0.1d0) then
print *, ' NOTE: Bi < 0.1'
print *, ' The lumped system analysis could also be'
print *, ' used for this problem (simpler method).'
print *, ''
end if
! =============================================
! FIND FIRST EIGENVALUE
! =============================================
call find_eigenvalue(geom_type, biot, zeta1)
call compute_coefficient(geom_type, zeta1, coeff_c1)
print *, '================================================='
print *, ' EIGENVALUE ANALYSIS'
print *, '================================================='
print *, ''
select case (geom_type)
case (1)
print *, ' Transcendental Equation:'
print *, ' zeta * tan(zeta) = Bi'
case (2)
print *, ' Transcendental Equation:'
print *, ' zeta * J1(zeta) / J0(zeta) = Bi'
case (3)
print *, ' Transcendental Equation:'
print *, ' 1 - zeta * cot(zeta) = Bi'
end select
print *, ''
print '(A,F14.8)', ' First eigenvalue (zeta1): ', zeta1
print '(A,F14.8)', ' Coefficient (C1): ', coeff_c1
print *, ''
! =============================================
! TEMPERATURE AT SPECIFIED POINT AND TIME
! =============================================
if (time_val > 0.0d0) then
! Center temperature
theta_center = coeff_c1 * exp(-(zeta1**2) * fourier)
! Temperature at position
call compute_theta_position(geom_type, zeta1, coeff_c1, fourier, x_ratio, theta_star)
temp_result = temp_inf + (temp_init - temp_inf) * theta_star
print *, '================================================='
print *, ' TEMPERATURE RESULTS'
print *, '================================================='
print *, ''
print *, ' One-term approximation:'
print *, ' theta* = C1 * exp(-zeta1^2 * Fo) * f(x)'
print *, ''
print '(A,F12.6)', ' theta* (center): ', theta_center
print '(A,F12.6)', ' theta* (at x/L or r/r0): ', theta_star
print '(A,F12.2,A)', ' T (center): ', &
temp_inf + (temp_init - temp_inf) * theta_center, ' deg-C'
print '(A,F12.2,A)', ' T (at position): ', temp_result, ' deg-C'
print '(A,F12.2,A)', ' T (surface, x=1): ', &
temp_inf + (temp_init - temp_inf) * theta_center * &
spatial_factor(geom_type, zeta1, 1.0d0), ' deg-C'
print *, ''
! Energy transfer
call compute_energy_ratio(geom_type, zeta1, coeff_c1, fourier, q_ratio)
print *, '================================================='
print *, ' ENERGY TRANSFER'
print *, '================================================='
print *, ''
print '(A,F12.6)', ' Q/Qmax: ', q_ratio
print '(A,F12.2,A)', ' Percentage: ', q_ratio * 100.0d0, ' %'
print *, ''
print *, ' Q/Qmax represents the fraction of maximum'
print *, ' possible energy transfer that has occurred'
print *, ' by time t.'
print *, ''
end if
! =============================================
! TEMPERATURE PROFILE AT TIME t
! =============================================
if (time_val > 0.0d0) then
print *, '================================================='
print *, ' TEMPERATURE PROFILE AT t =', time_val, ' s'
print *, '================================================='
print *, ''
select case (geom_type)
case (1)
print *, ' x/L | T (deg-C) | theta*'
case default
print *, ' r/r0 | T (deg-C) | theta*'
end select
print *, ' -------------------------------------------'
n_points = 20
do i = 0, n_points
x_pos = dble(i) / dble(n_points)
call compute_theta_position(geom_type, zeta1, coeff_c1, fourier, x_pos, theta_star)
temp_x = temp_inf + (temp_init - temp_inf) * theta_star
write(*, '(3X,F8.4,5X,A,2X,F10.3,5X,A,2X,F10.6)') &
x_pos, '|', temp_x, '|', theta_star
end do
print *, ''
end if
! =============================================
! TEMPERATURE EVOLUTION AT CENTER
! =============================================
print *, '================================================='
print *, ' TEMPERATURE EVOLUTION AT CENTER'
print *, '================================================='
print *, ''
print *, ' Fo | t (s) | T_center | theta*'
print *, ' --------------------------------------------------------'
n_points = 40
delta_fo = 5.0d0 / dble(n_points)
do i = 0, n_points
fo_val = dble(i) * delta_fo
if (fo_val < 0.001d0) fo_val = 0.001d0
theta_center = coeff_c1 * exp(-(zeta1**2) * fo_val)
temp_x = temp_inf + (temp_init - temp_inf) * theta_center
write(*, '(3X,F8.4,5X,A,2X,F12.2,3X,A,2X,F10.3,3X,A,2X,F10.6)') &
fo_val, '|', fo_val * dim_char**2 / alpha_diff, &
'|', temp_x, '|', theta_center
end do
print *, ''
! =============================================
! SURFACE TEMPERATURE EVOLUTION
! =============================================
print *, '================================================='
print *, ' SURFACE TEMPERATURE EVOLUTION'
print *, '================================================='
print *, ''
print *, ' Fo | t (s) | T_surface | theta*_s'
print *, ' --------------------------------------------------------'
do i = 0, n_points
fo_val = dble(i) * delta_fo
if (fo_val < 0.001d0) fo_val = 0.001d0
theta_center = coeff_c1 * exp(-(zeta1**2) * fo_val)
theta_star = theta_center * spatial_factor(geom_type, zeta1, 1.0d0)
temp_x = temp_inf + (temp_init - temp_inf) * theta_star
write(*, '(3X,F8.4,5X,A,2X,F12.2,3X,A,2X,F10.3,3X,A,2X,F10.6)') &
fo_val, '|', fo_val * dim_char**2 / alpha_diff, &
'|', temp_x, '|', theta_star
end do
print *, ''
! =============================================
! RECOMMENDATIONS
! =============================================
print *, '================================================='
print *, ' DESIGN RECOMMENDATIONS'
print *, '================================================='
print *, ''
if (biot < 0.1d0) then
print *, ' LOW Bi (< 0.1): Internal resistance negligible'
print *, ' - Lumped system analysis is also applicable'
print *, ' - Temperature is nearly uniform inside'
else if (biot < 1.0d0) then
print *, ' MODERATE Bi (0.1 - 1.0):'
print *, ' - Both internal and external resistances matter'
print *, ' - One-term approximation works well for Fo > 0.2'
else if (biot < 10.0d0) then
print *, ' HIGH Bi (1 - 10):'
print *, ' - Internal resistance dominates'
print *, ' - Significant temperature gradients inside'
else
print *, ' VERY HIGH Bi (> 10):'
print *, ' - Surface temperature approaches T_inf quickly'
print *, ' - Large internal temperature gradients'
print *, ' - Consider if surface boundary approaches'
print *, ' a constant temperature condition'
end if
print *, ''
print *, '================================================='
print *, ' CALCULATION COMPLETE'
print *, '================================================='
print *, ''
contains
! =========================================================
! BESSEL FUNCTION J0 (Taylor series, 20 terms)
! =========================================================
function bessel_j0(x) result(j0)
implicit none
real(8), intent(in) :: x
real(8) :: j0, term, x2
integer :: m
j0 = 1.0d0
term = 1.0d0
x2 = (x / 2.0d0)**2
do m = 1, 25
term = -term * x2 / (dble(m)**2)
j0 = j0 + term
if (abs(term) < 1.0d-15) exit
end do
end function bessel_j0
! =========================================================
! BESSEL FUNCTION J1 (Taylor series, 20 terms)
! =========================================================
function bessel_j1(x) result(j1)
implicit none
real(8), intent(in) :: x
real(8) :: j1, term, x2
integer :: m
j1 = x / 2.0d0
term = x / 2.0d0
x2 = (x / 2.0d0)**2
do m = 1, 25
term = -term * x2 / (dble(m) * dble(m + 1))
j1 = j1 + term
if (abs(term) < 1.0d-15) exit
end do
end function bessel_j1
! =========================================================
! SPATIAL FACTOR at position x_r (x/L or r/r0)
! =========================================================
function spatial_factor(gtype, zeta, x_r) result(sf)
implicit none
integer, intent(in) :: gtype
real(8), intent(in) :: zeta, x_r
real(8) :: sf
select case (gtype)
case (1) ! Plane wall: cos(zeta * x/L)
sf = cos(zeta * x_r)
case (2) ! Cylinder: J0(zeta * r/r0)
sf = bessel_j0(zeta * x_r)
case (3) ! Sphere: sin(zeta * r/r0) / (zeta * r/r0)
if (x_r < 1.0d-10) then
sf = 1.0d0 ! limit as r->0
else
sf = sin(zeta * x_r) / (zeta * x_r)
end if
end select
end function spatial_factor
! =========================================================
! FIND FIRST EIGENVALUE via Newton-Raphson
! =========================================================
subroutine find_eigenvalue(gtype, bi, zeta)
implicit none
integer, intent(in) :: gtype
real(8), intent(in) :: bi
real(8), intent(out) :: zeta
real(8) :: f_val, df_val, correction
real(8) :: j0v, j1v
integer :: iter
! Initial guess
select case (gtype)
case (1) ! zeta * tan(zeta) = Bi
if (bi < 0.01d0) then
zeta = sqrt(bi)
else if (bi < 1.0d0) then
zeta = sqrt(bi) * (1.0d0 - bi/6.0d0)
else if (bi < 10.0d0) then
zeta = 1.0d0 + 0.07d0 * bi
else
zeta = pi/2.0d0 - 0.1d0
end if
case (2) ! zeta * J1(zeta) / J0(zeta) = Bi
if (bi < 0.01d0) then
zeta = sqrt(2.0d0 * bi)
else if (bi < 1.0d0) then
zeta = sqrt(2.0d0 * bi)
else if (bi < 10.0d0) then
zeta = 1.5d0 + 0.05d0 * bi
else
zeta = 2.4048d0 - 0.1d0
end if
case (3) ! 1 - zeta * cot(zeta) = Bi
if (bi < 0.01d0) then
zeta = sqrt(3.0d0 * bi)
else if (bi < 1.0d0) then
zeta = sqrt(3.0d0 * bi)
else if (bi < 10.0d0) then
zeta = 2.0d0 + 0.05d0 * bi
else
zeta = pi - 0.1d0
end if
end select
! Newton-Raphson iterations
do iter = 1, 200
select case (gtype)
case (1) ! f = zeta * tan(zeta) - Bi
if (abs(cos(zeta)) < 1.0d-12) then
zeta = zeta - 0.01d0
cycle
end if
f_val = zeta * tan(zeta) - bi
df_val = tan(zeta) + zeta / (cos(zeta)**2)
case (2) ! f = zeta * J1(zeta) / J0(zeta) - Bi
j0v = bessel_j0(zeta)
j1v = bessel_j1(zeta)
if (abs(j0v) < 1.0d-12) then
zeta = zeta - 0.01d0
cycle
end if
f_val = zeta * j1v / j0v - bi
! Derivative using Bessel recurrences
df_val = j1v/j0v + zeta * (j0v * (j0v/zeta - j1v) - &
j1v * (-j1v)) / (j0v**2)
case (3) ! f = 1 - zeta * cos(zeta)/sin(zeta) - Bi
if (abs(sin(zeta)) < 1.0d-12) then
zeta = zeta - 0.01d0
cycle
end if
f_val = 1.0d0 - zeta * cos(zeta) / sin(zeta) - bi
df_val = -cos(zeta)/sin(zeta) + zeta/(sin(zeta)**2)
end select
if (abs(df_val) < 1.0d-15) exit
correction = f_val / df_val
! Damping for stability
if (abs(correction) > 0.5d0) then
correction = sign(0.5d0, correction)
end if
zeta = zeta - correction
! Keep zeta positive and in first root range
if (zeta <= 0.0d0) zeta = 0.01d0
select case (gtype)
case (1)
if (zeta >= pi/2.0d0) zeta = pi/2.0d0 - 0.001d0
case (2)
if (zeta >= 2.4048d0) zeta = 2.4048d0 - 0.001d0
case (3)
if (zeta >= pi) zeta = pi - 0.001d0
end select
if (abs(correction) < 1.0d-12) exit
end do
end subroutine find_eigenvalue
! =========================================================
! COMPUTE C1 COEFFICIENT
! =========================================================
subroutine compute_coefficient(gtype, zeta, c1)
implicit none
integer, intent(in) :: gtype
real(8), intent(in) :: zeta
real(8), intent(out) :: c1
select case (gtype)
case (1) ! C1 = 4*sin(zeta) / (2*zeta + sin(2*zeta))
c1 = 4.0d0 * sin(zeta) / (2.0d0 * zeta + sin(2.0d0 * zeta))
case (2) ! C1 = 2 * J1(zeta) / (zeta * (J0^2 + J1^2))
c1 = 2.0d0 / zeta * bessel_j1(zeta) / &
(bessel_j0(zeta)**2 + bessel_j1(zeta)**2)
case (3) ! C1 = 4*(sin(zeta) - zeta*cos(zeta)) / (2*zeta - sin(2*zeta))
c1 = 4.0d0 * (sin(zeta) - zeta * cos(zeta)) / &
(2.0d0 * zeta - sin(2.0d0 * zeta))
end select
end subroutine compute_coefficient
! =========================================================
! COMPUTE THETA* AT A POSITION
! =========================================================
subroutine compute_theta_position(gtype, zeta, c1, fo, x_r, theta)
implicit none
integer, intent(in) :: gtype
real(8), intent(in) :: zeta, c1, fo, x_r
real(8), intent(out) :: theta
theta = c1 * exp(-(zeta**2) * fo) * spatial_factor(gtype, zeta, x_r)
! Clamp to valid range
if (theta < 0.0d0) theta = 0.0d0
if (theta > 1.0d0) theta = 1.0d0
end subroutine compute_theta_position
! =========================================================
! COMPUTE ENERGY RATIO Q/Qmax
! =========================================================
subroutine compute_energy_ratio(gtype, zeta, c1, fo, qr)
implicit none
integer, intent(in) :: gtype
real(8), intent(in) :: zeta, c1, fo
real(8), intent(out) :: qr
real(8) :: theta0
theta0 = c1 * exp(-(zeta**2) * fo)
select case (gtype)
case (1) ! Q/Qmax = 1 - theta0*sin(zeta)/zeta
qr = 1.0d0 - theta0 * sin(zeta) / zeta
case (2) ! Q/Qmax = 1 - 2*theta0*J1(zeta)/zeta
qr = 1.0d0 - 2.0d0 * theta0 * bessel_j1(zeta) / zeta
case (3) ! Q/Qmax = 1 - 3*theta0*(sin(zeta)-zeta*cos(zeta))/zeta^3
qr = 1.0d0 - 3.0d0 * theta0 * &
(sin(zeta) - zeta * cos(zeta)) / (zeta**3)
end select
if (qr < 0.0d0) qr = 0.0d0
if (qr > 1.0d0) qr = 1.0d0
end subroutine compute_energy_ratio
end program transient_conduction_1d
Solver Description
Calculates transient temperature distribution in plane walls, cylinders, or spheres using the one-term approximation of Fourier series. Solves transcendental boundary conditions to compute first eigenvalues ($\lambda_1$) and coefficients ($A_1$) dynamically. Valid for Fourier numbers $Fo > 0.2$.
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 transient_conduction_1d.f90 -o transient_1d
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
transient_1d < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Geometry Code (1=Plane wall, 2=Cylinder, 3=Sphere)
3
! Half-Thickness or Outer Radius ($L$) [m]
0.06
! Conductivity ($k$) [W/m-K]
110.0
! Density ($\rho$) [kg/m³]
8530.0
! Specific Heat ($c_p$) [J/kg-K]
380.0
! Convection Coefficient ($h$) [W/m²-K]
120.0
! Initial Temperature ($T_i$) [°C]
150.0
! Ambient Temperature ($T_\infty$) [°C]
20.0
! Normalized Position ($x/L$ or $r/r_0$) [0-1]
0.0
! Time ($t$) [s]
2700.0
3
! Half-Thickness or Outer Radius ($L$) [m]
0.06
! Conductivity ($k$) [W/m-K]
110.0
! Density ($\rho$) [kg/m³]
8530.0
! Specific Heat ($c_p$) [J/kg-K]
380.0
! Convection Coefficient ($h$) [W/m²-K]
120.0
! Initial Temperature ($T_i$) [°C]
150.0
! Ambient Temperature ($T_\infty$) [°C]
20.0
! Normalized Position ($x/L$ or $r/r_0$) [0-1]
0.0
! Time ($t$) [s]
2700.0