💻 Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Transient Conduction — 1D Series
Core Numerical Engine in Fortran 90 • 33 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
For transient conduction in 1D shapes (plane walls, long cylinders, spheres) where the Biot number is greater than 0.1, internal temperature gradients are non-negligible. The one-term Heisler approximation is 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 -O3 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 (1=Slab, 2=Cylinder, 3=Sphere)
1
! Characteristic length/radius [m]
0.06
! Thermal conductivity k [W/m-K]
237.0
! Density rho [kg/m3]
2702.0
! Specific heat Cp [J/kg-K]
903.0
! Convection coefficient h [W/m2-K]
500.0
! Initial temperature Ti [°C]
200.0
! Ambient temperature Tinf [°C]
40.0
! Position ratio x/L or r/ro
0.0
! Time t [s]
60.0
1
! Characteristic length/radius [m]
0.06
! Thermal conductivity k [W/m-K]
237.0
! Density rho [kg/m3]
2702.0
! Specific heat Cp [J/kg-K]
903.0
! Convection coefficient h [W/m2-K]
500.0
! Initial temperature Ti [°C]
200.0
! Ambient temperature Tinf [°C]
40.0
! Position ratio x/L or r/ro
0.0
! Time t [s]
60.0