💻 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.
Multidimensional Transient Conduction
Core Numerical Engine in Fortran 90 • 45 total downloads
multidimensional_transient.f90
! =========================================================================
! Source File: multidimensional_transient.f90
! =========================================================================
program multidimensional_transient
implicit none
! Variable declarations
integer :: geom_type, i
real(8) :: k_mat, rho_mat, cp_mat, alpha_diff
real(8) :: h1, h2, h3
real(8) :: temp_init, temp_inf
real(8) :: dim1, dim2, dim3
real(8) :: pos1, pos2, pos3
real(8) :: time_val
! 1D results for components
real(8) :: bi1, bi2, bi3
real(8) :: fo1, fo2, fo3
real(8) :: zeta1_1, zeta1_2, zeta1_3
real(8) :: c1_1, c1_2, c1_3
real(8) :: theta1, theta2, theta3
real(8) :: q_ratio1, q_ratio2, q_ratio3
! Combined results
real(8) :: theta_combined, temp_combined
real(8) :: q_ratio_combined
real(8), parameter :: pi = 3.141592653589793d0
! Read inputs
read(*,*) geom_type ! 1=Short Cylinder, 2=Rectangular Bar, 3=Rectangular Block, 4=Semi-infinite Cylinder, 5=Semi-infinite Rectangular Bar
read(*,*) k_mat ! Thermal conductivity [W/m.K]
read(*,*) rho_mat ! Density [kg/m3]
read(*,*) cp_mat ! Specific heat [J/kg.K]
read(*,*) h1, h2, h3 ! Convection coefficients for directions 1, 2, 3
read(*,*) temp_init ! Initial temperature [deg-C]
read(*,*) temp_inf ! Ambient fluid temperature [deg-C]
read(*,*) dim1, dim2, dim3 ! Specific geometry dimensions (m)
read(*,*) pos1, pos2, pos3 ! Specific coordinates (m)
read(*,*) time_val ! Time elapsed [s]
! Safety limits
if (k_mat <= 0.0d0) k_mat = 1.0d0
if (rho_mat <= 0.0d0) rho_mat = 1000.0d0
if (cp_mat <= 0.0d0) cp_mat = 1000.0d0
if (time_val <= 0.0d0) time_val = 1.0d-5
! Thermal diffusivity
alpha_diff = k_mat / (rho_mat * cp_mat)
! Initialize 1D component variables
bi1 = 0.0d0; bi2 = 0.0d0; bi3 = 0.0d0
fo1 = 0.0d0; fo2 = 0.0d0; fo3 = 0.0d0
zeta1_1 = 0.0d0; zeta1_2 = 0.0d0; zeta1_3 = 0.0d0
c1_1 = 0.0d0; c1_2 = 0.0d0; c1_3 = 0.0d0
theta1 = 1.0d0; theta2 = 1.0d0; theta3 = 1.0d0
q_ratio1 = 0.0d0; q_ratio2 = 0.0d0; q_ratio3 = 0.0d0
! Compute 1D components based on combined geometry choice
select case (geom_type)
case (1) ! Short Cylinder = Plane Wall (dim1) x Infinite Cylinder (dim2)
! Component 1: Plane Wall
if (dim1 > 0.0d0) then
bi1 = h1 * dim1 / k_mat
fo1 = alpha_diff * time_val / (dim1**2)
call find_eigenvalue(1, bi1, zeta1_1)
call compute_coefficient(1, zeta1_1, c1_1)
call compute_theta_position(1, zeta1_1, c1_1, fo1, pos1/dim1, theta1)
call compute_energy_ratio(1, zeta1_1, c1_1, fo1, q_ratio1)
end if
! Component 2: Infinite Cylinder
if (dim2 > 0.0d0) then
bi2 = h2 * dim2 / k_mat
fo2 = alpha_diff * time_val / (dim2**2)
call find_eigenvalue(2, bi2, zeta1_2)
call compute_coefficient(2, zeta1_2, c1_2)
call compute_theta_position(2, zeta1_2, c1_2, fo2, pos2/dim2, theta2)
call compute_energy_ratio(2, zeta1_2, c1_2, fo2, q_ratio2)
end if
! Combine
theta_combined = theta1 * theta2
q_ratio_combined = 1.0d0 - (1.0d0 - q_ratio1) * (1.0d0 - q_ratio2)
case (2) ! Rectangular Bar = Wall 1 (dim1) x Wall 2 (dim2)
! Component 1: Wall 1
if (dim1 > 0.0d0) then
bi1 = h1 * dim1 / k_mat
fo1 = alpha_diff * time_val / (dim1**2)
call find_eigenvalue(1, bi1, zeta1_1)
call compute_coefficient(1, zeta1_1, c1_1)
call compute_theta_position(1, zeta1_1, c1_1, fo1, pos1/dim1, theta1)
call compute_energy_ratio(1, zeta1_1, c1_1, fo1, q_ratio1)
end if
! Component 2: Wall 2
if (dim2 > 0.0d0) then
bi2 = h2 * dim2 / k_mat
fo2 = alpha_diff * time_val / (dim2**2)
call find_eigenvalue(1, bi2, zeta1_2)
call compute_coefficient(1, zeta1_2, c1_2)
call compute_theta_position(1, zeta1_2, c1_2, fo2, pos2/dim2, theta2)
call compute_energy_ratio(1, zeta1_2, c1_2, fo2, q_ratio2)
end if
! Combine
theta_combined = theta1 * theta2
q_ratio_combined = 1.0d0 - (1.0d0 - q_ratio1) * (1.0d0 - q_ratio2)
case (3) ! Rectangular Block = Wall 1 (dim1) x Wall 2 (dim2) x Wall 3 (dim3)
! Component 1: Wall 1
if (dim1 > 0.0d0) then
bi1 = h1 * dim1 / k_mat
fo1 = alpha_diff * time_val / (dim1**2)
call find_eigenvalue(1, bi1, zeta1_1)
call compute_coefficient(1, zeta1_1, c1_1)
call compute_theta_position(1, zeta1_1, c1_1, fo1, pos1/dim1, theta1)
call compute_energy_ratio(1, zeta1_1, c1_1, fo1, q_ratio1)
end if
! Component 2: Wall 2
if (dim2 > 0.0d0) then
bi2 = h2 * dim2 / k_mat
fo2 = alpha_diff * time_val / (dim2**2)
call find_eigenvalue(1, bi2, zeta1_2)
call compute_coefficient(1, zeta1_2, c1_2)
call compute_theta_position(1, zeta1_2, c1_2, fo2, pos2/dim2, theta2)
call compute_energy_ratio(1, zeta1_2, c1_2, fo2, q_ratio2)
end if
! Component 3: Wall 3
if (dim3 > 0.0d0) then
bi3 = h3 * dim3 / k_mat
fo3 = alpha_diff * time_val / (dim3**2)
call find_eigenvalue(1, bi3, zeta1_3)
call compute_coefficient(1, zeta1_3, c1_3)
call compute_theta_position(1, zeta1_3, c1_3, fo3, pos3/dim3, theta3)
call compute_energy_ratio(1, zeta1_3, c1_3, fo3, q_ratio3)
end if
! Combine
theta_combined = theta1 * theta2 * theta3
q_ratio_combined = 1.0d0 - (1.0d0 - q_ratio1) * (1.0d0 - q_ratio2) * (1.0d0 - q_ratio3)
case (4) ! Semi-infinite Cylinder = Semi-infinite Solid (pos1) x Infinite Cylinder (dim2)
! Component 1: Semi-infinite Solid (evaluates at depth pos1)
call compute_semi_infinite(pos1, h1, k_mat, alpha_diff, time_val, theta1)
q_ratio1 = -1.0d0 ! Not applicable
! Component 2: Infinite Cylinder (dim2)
if (dim2 > 0.0d0) then
bi2 = h2 * dim2 / k_mat
fo2 = alpha_diff * time_val / (dim2**2)
call find_eigenvalue(2, bi2, zeta1_2)
call compute_coefficient(2, zeta1_2, c1_2)
call compute_theta_position(2, zeta1_2, c1_2, fo2, pos2/dim2, theta2)
call compute_energy_ratio(2, zeta1_2, c1_2, fo2, q_ratio2)
end if
! Combine
theta_combined = theta1 * theta2
q_ratio_combined = -1.0d0 ! Infinite Volume: N/A
case (5) ! Semi-infinite Rectangular Bar = Semi-infinite Solid (pos1) x Plane Wall (dim2)
! Component 1: Semi-infinite Solid (evaluates at depth pos1)
call compute_semi_infinite(pos1, h1, k_mat, alpha_diff, time_val, theta1)
q_ratio1 = -1.0d0 ! Not applicable
! Component 2: Plane Wall (dim2)
if (dim2 > 0.0d0) then
bi2 = h2 * dim2 / k_mat
fo2 = alpha_diff * time_val / (dim2**2)
call find_eigenvalue(1, bi2, zeta1_2)
call compute_coefficient(1, zeta1_2, c1_2)
call compute_theta_position(1, zeta1_2, c1_2, fo2, pos2/dim2, theta2)
call compute_energy_ratio(1, zeta1_2, c1_2, fo2, q_ratio2)
end if
! Combine
theta_combined = theta1 * theta2
q_ratio_combined = -1.0d0 ! Infinite Volume: N/A
case default
print *, "ERROR: Invalid geometry type selection"
stop
end select
! Resulting Temperature
temp_combined = temp_inf + (temp_init - temp_inf) * theta_combined
! ==================================================
! DISPLAY REPORT
! ==================================================
print *, "=================================================="
print *, " MULTIDIMENSIONAL TRANSIENT CONDUCTION REPORT"
print *, " (Product Solution)"
print *, "=================================================="
print *, ""
print *, "1. THERMOPHYSICAL PROPERTIES"
print *, "--------------------------------------------------"
print '(A, F10.4, A)', " Conductivity (k): ", k_mat, " W/m.K"
print '(A, F10.2, A)', " Density (rho): ", rho_mat, " kg/m3"
print '(A, F10.2, A)', " Specific Heat (cp): ", cp_mat, " J/kg.K"
print '(A, ES12.4, A)'," Thermal Diffusivity (a): ", alpha_diff, " m2/s"
print '(A, F10.2, A)', " Initial Temp (Ti): ", temp_init, " C"
print '(A, F10.2, A)', " Ambient Fluid Temp (Tinf):", temp_inf, " C"
print *, ""
print *, "2. GEOMETRY CONFIGURATION"
print *, "--------------------------------------------------"
select case (geom_type)
case (1)
print *, " Type: Short Cylinder (Plane Wall x Infinite Cylinder)"
print '(A, F10.4, A)', " Wall Half-Thickness (L1):", dim1, " m"
print '(A, F10.4, A)', " Cylinder Radius (r0): ", dim2, " m"
print '(A, F10.4, A)', " Axial Coordinate (x): ", pos1, " m"
print '(A, F10.4, A)', " Radial Coordinate (r): ", pos2, " m"
print '(A, F10.2, A)', " Wall Convection (h1): ", h1, " W/m2.K"
print '(A, F10.2, A)', " Cylinder Convection (h2):", h2, " W/m2.K"
case (2)
print *, " Type: Rectangular Bar (Wall 1 x Wall 2)"
print '(A, F10.4, A)', " Wall 1 Half-Thickness(L1):", dim1, " m"
print '(A, F10.4, A)', " Wall 2 Half-Thickness(L2):", dim2, " m"
print '(A, F10.4, A)', " Coordinate x1: ", pos1, " m"
print '(A, F10.4, A)', " Coordinate x2: ", pos2, " m"
print '(A, F10.2, A)', " Convection Face 1 (h1): ", h1, " W/m2.K"
print '(A, F10.2, A)', " Convection Face 2 (h2): ", h2, " W/m2.K"
case (3)
print *, " Type: Rectangular Block (Wall 1 x Wall 2 x Wall 3)"
print '(A, F10.4, A)', " Wall 1 Half-Thickness(L1):", dim1, " m"
print '(A, F10.4, A)', " Wall 2 Half-Thickness(L2):", dim2, " m"
print '(A, F10.4, A)', " Wall 3 Half-Thickness(L3):", dim3, " m"
print '(A, F10.4, A)', " Coordinate x1: ", pos1, " m"
print '(A, F10.4, A)', " Coordinate x2: ", pos2, " m"
print '(A, F10.4, A)', " Coordinate x3: ", pos3, " m"
print '(A, F10.2, A)', " Convection Face 1 (h1): ", h1, " W/m2.K"
print '(A, F10.2, A)', " Convection Face 2 (h2): ", h2, " W/m2.K"
print '(A, F10.2, A)', " Convection Face 3 (h3): ", h3, " W/m2.K"
case (4)
print *, " Type: Semi-Infinite Cylinder (Semi-Infinite Solid x Infinite Cylinder)"
print '(A, F10.4, A)', " Cylinder Radius (r0): ", dim2, " m"
print '(A, F10.4, A)', " Semi-Infinite Depth (x): ", pos1, " m"
print '(A, F10.4, A)', " Radial Coordinate (r): ", pos2, " m"
print '(A, F10.2, A)', " End Convection (h1): ", h1, " W/m2.K"
print '(A, F10.2, A)', " Cylinder Convection (h2):", h2, " W/m2.K"
case (5)
print *, " Type: Semi-Infinite Rectangular Bar (Semi-Infinite Solid x Plane Wall)"
print '(A, F10.4, A)', " Wall Half-Thickness (L1):", dim2, " m"
print '(A, F10.4, A)', " Semi-Infinite Depth (x): ", pos1, " m"
print '(A, F10.4, A)', " Wall Coordinate (y): ", pos2, " m"
print '(A, F10.2, A)', " End Convection (h1): ", h1, " W/m2.K"
print '(A, F10.2, A)', " Wall Convection (h2): ", h2, " W/m2.K"
end select
print '(A, F10.2, A)', " Elapsed Time (t): ", time_val, " s"
print *, ""
print *, "3. 1D SUB-PROBLEM ANALYSIS"
print *, "--------------------------------------------------"
select case (geom_type)
case (1)
print *, " [Sub-Problem 1: Plane Wall]"
call print_1d_metrics(bi1, fo1, zeta1_1, c1_1, theta1)
print *, ""
print *, " [Sub-Problem 2: Infinite Cylinder]"
call print_1d_metrics(bi2, fo2, zeta1_2, c1_2, theta2)
case (2)
print *, " [Sub-Problem 1: Plane Wall 1]"
call print_1d_metrics(bi1, fo1, zeta1_1, c1_1, theta1)
print *, ""
print *, " [Sub-Problem 2: Plane Wall 2]"
call print_1d_metrics(bi2, fo2, zeta1_2, c1_2, theta2)
case (3)
print *, " [Sub-Problem 1: Plane Wall 1]"
call print_1d_metrics(bi1, fo1, zeta1_1, c1_1, theta1)
print *, ""
print *, " [Sub-Problem 2: Plane Wall 2]"
call print_1d_metrics(bi2, fo2, zeta1_2, c1_2, theta2)
print *, ""
print *, " [Sub-Problem 3: Plane Wall 3]"
call print_1d_metrics(bi3, fo3, zeta1_3, c1_3, theta3)
case (4)
print *, " [Sub-Problem 1: Semi-Infinite Solid]"
print '(A, F12.6)', " Dimensionless Temp (theta1): ", theta1
print *, ""
print *, " [Sub-Problem 2: Infinite Cylinder]"
call print_1d_metrics(bi2, fo2, zeta1_2, c1_2, theta2)
case (5)
print *, " [Sub-Problem 1: Semi-Infinite Solid]"
print '(A, F12.6)', " Dimensionless Temp (theta1): ", theta1
print *, ""
print *, " [Sub-Problem 2: Plane Wall]"
call print_1d_metrics(bi2, fo2, zeta1_2, c1_2, theta2)
end select
print *, ""
print *, "4. COMBINED RESULTS (PRODUCT SOLUTION)"
print *, "--------------------------------------------------"
print '(A, F12.6)', " Dimensionless Temp (theta*): ", theta_combined
print '(A, F12.2, A)', " Combined Temperature T(pos, t): ", temp_combined, " C"
if (q_ratio_combined >= 0.0d0) then
print '(A, F12.6)', " Heat Transfer Ratio Q/Qmax: ", q_ratio_combined
print '(A, F12.2, A)', " Percentage Q/Qmax: ", q_ratio_combined * 100.0d0, " %"
else
print *, " Heat Transfer Ratio Q/Qmax: N/A (Infinite Volume)"
end if
print *, ""
print *, "=================================================="
print *, " CALCULATION COMPLETE"
print *, "=================================================="
contains
! =========================================================
! PRINTS 1D SUB-PROBLEM METRICS WITH ACCURACY CHECKS
! =========================================================
subroutine print_1d_metrics(bi, fo, zeta, c, theta)
real(8), intent(in) :: bi, fo, zeta, c, theta
print '(A, F12.4)', " Biot Number (Bi): ", bi
print '(A, F12.4)', " Fourier Number (Fo): ", fo
print '(A, F12.6)', " First Eigenvalue (zeta1): ", zeta
print '(A, F12.6)', " Coefficient (C1): ", c
print '(A, F12.6)', " Dimensionless Temp (theta): ", theta
if (fo < 0.2d0) then
print *, " WARNING: Fo < 0.2. One-term approximation may be inaccurate."
else
print *, " STATUS: Fo >= 0.2. One-term approximation is valid."
end if
end subroutine print_1d_metrics
! =========================================================
! STABLE EXP(U^2) * ERFC(U) TO AVOID FLOATING OVERFLOW
! =========================================================
function exp_erfc_stable(u) result(val)
real(8), intent(in) :: u
real(8) :: val, u2
if (u > 10.0d0) then
u2 = u * u
! Asymptotic series expansion
val = 1.0d0 / (u * sqrt(pi)) * (1.0d0 - 0.5d0 / u2 + 0.75d0 / (u2 * u2) - 1.875d0 / (u2 * u2 * u2))
else
val = exp(u * u) * erfc(u)
end if
end function exp_erfc_stable
! =========================================================
! SOLVES THE SEMI-INFINITE CONVECTION SUB-PROBLEM
! =========================================================
subroutine compute_semi_infinite(x, h, k, alpha, t, theta)
real(8), intent(in) :: x, h, k, alpha, t
real(8), intent(out) :: theta
real(8) :: arg1, arg2
arg1 = x / (2.0d0 * sqrt(alpha * t))
arg2 = h * sqrt(alpha * t) / k
! Math equation: theta = erf(arg1) + exp(h*x/k + h^2*alpha*t/k^2) * erfc(arg1 + arg2)
! Stably computed as: erf(arg1) + exp(-arg1^2) * exp_erfc_stable(arg1 + arg2)
theta = erf(arg1) + exp(-arg1**2) * exp_erfc_stable(arg1 + arg2)
if (theta < 0.0d0) theta = 0.0d0
if (theta > 1.0d0) theta = 1.0d0
end subroutine compute_semi_infinite
! =========================================================
! BESSEL FUNCTION J0 (Taylor series, 25 terms)
! =========================================================
function bessel_j0(x) result(j0)
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, 25 terms)
! =========================================================
function bessel_j1(x) result(j1)
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)
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)
end select
end function spatial_factor
! =========================================================
! FIND FIRST EIGENVALUE via Newton-Raphson
! =========================================================
subroutine find_eigenvalue(gtype, bi, zeta)
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
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)
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
end select
if (abs(correction) < 1.0d-12) exit
end do
end subroutine find_eigenvalue
! =========================================================
! COMPUTE C1 COEFFICIENT
! =========================================================
subroutine compute_coefficient(gtype, zeta, c1)
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)
end select
end subroutine compute_coefficient
! =========================================================
! COMPUTE THETA* AT A POSITION
! =========================================================
subroutine compute_theta_position(gtype, zeta, c1, fo, x_r, theta)
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)
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
end select
if (qr < 0.0d0) qr = 0.0d0
if (qr > 1.0d0) qr = 1.0d0
end subroutine compute_energy_ratio
end program multidimensional_transient
Solver Description
Solve multi-dimensional transient heat conduction problems using the product solution method. Interactive 2D cross-sectional temperature profile heatmaps and detailed Fortran calculations.
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 multidimensional_transient.f90 -o multidimensional_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
multidimensional_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Geometry type (1=Rectangular block, 2=Cylinder block)
1
! Thermal conductivity k [W/m-K]
43.0
! Density rho [kg/m3]
7850.0
! Specific heat Cp [J/kg-K]
475.0
! Convection coefficients h1 h2 h3 [W/m2-K]
120.0 200.0 150.0
! Initial temperature Ti [°C]
20.0
! Ambient fluid temperature T_inf [°C]
180.0
! Half-thickness dimensions x1 x2 x3 [m]
0.05 0.03 0.04
! Evaluation coordinates x y z [m]
0.0 0.0 0.0
! Time t [s]
120.0
1
! Thermal conductivity k [W/m-K]
43.0
! Density rho [kg/m3]
7850.0
! Specific heat Cp [J/kg-K]
475.0
! Convection coefficients h1 h2 h3 [W/m2-K]
120.0 200.0 150.0
! Initial temperature Ti [°C]
20.0
! Ambient fluid temperature T_inf [°C]
180.0
! Half-thickness dimensions x1 x2 x3 [m]
0.05 0.03 0.04
! Evaluation coordinates x y z [m]
0.0 0.0 0.0
! Time t [s]
120.0