π¬
Solver Purpose & Physical Scope
Calculate 2D multidimensional transient heat conduction in rectangular billets: center core temperature (T0), corner surface temperature, Biot & Fourier numbers, and heat removed (Q/Qmax).
π Discipline: Conduction
β‘ Precision: IEEE-754 64-bit Real(`real(8)`)
π₯ Total Downloads: 321 times
π Source File:
transient_multidim_conduction_adi.f90
π calcul/Conduction /
transient_multidim_conduction_adi.f90
program transient_multidim_conduction_adi
implicit none
integer :: iostat_val, iter
double precision :: Lx_mm, Ly_mm, k_cond, rho_dens, cp_spec, Ti_C, Tinf_C, h_conv, time_sec
double precision :: Lx_m, Ly_m, alpha_diff, Bix, Biy, Fox, Foy
double precision :: zeta_x, zeta_y, C_x, C_y, theta_x0, theta_y0, theta_x_surf, theta_y_surf
double precision :: theta_center, theta_corner, T_center_C, T_corner_C, Q_ratio
double precision :: z_low, z_high, z_mid, f_mid
double precision, parameter :: PI = 3.141592653589793d0
! Read inputs
read(*,*,iostat=iostat_val) Lx_mm ! Half-width in X [mm] (e.g. 50.0)
read(*,*,iostat=iostat_val) Ly_mm ! Half-width in Y [mm] (e.g. 80.0)
read(*,*,iostat=iostat_val) k_cond ! Thermal Conductivity [W/(m.K)] (e.g. 45.0 for carbon steel)
read(*,*,iostat=iostat_val) rho_dens ! Density [kg/m3] (e.g. 7850.0)
read(*,*,iostat=iostat_val) cp_spec ! Specific Heat [J/(kg.K)] (e.g. 480.0)
read(*,*,iostat=iostat_val) Ti_C ! Initial Temperature [deg C] (e.g. 850.0)
read(*,*,iostat=iostat_val) Tinf_C ! Ambient Quench Temp [deg C] (e.g. 30.0)
read(*,*,iostat=iostat_val) h_conv ! Convection Coeff [W/(m2.K)] (e.g. 350.0)
read(*,*,iostat=iostat_val) time_sec ! Elapsed Time [s] (e.g. 300.0)
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid input data for 2D transient conduction calculation.'
stop
end if
if (Lx_mm <= 0.0d0 .or. Ly_mm <= 0.0d0 .or. k_cond <= 0.0d0 .or. time_sec < 0.0d0) then
write(*,*) 'ERROR: Dimensions, conductivity, and time must be positive.'
stop
end if
Lx_m = Lx_mm * 1.0d-3
Ly_m = Ly_mm * 1.0d-3
alpha_diff = k_cond / (rho_dens * cp_spec) ! m2/s
Bix = (h_conv * Lx_m) / k_cond
Biy = (h_conv * Ly_m) / k_cond
Fox = (alpha_diff * time_sec) / (Lx_m**2)
Foy = (alpha_diff * time_sec) / (Ly_m**2)
! Solve transcendental root zeta_x * tan(zeta_x) = Bix
z_low = 0.0001d0
z_high = min(1.5707d0, max(0.1d0, atan(Bix)))
if (Bix > 10.0d0) z_high = 1.5700d0
do iter = 1, 60
z_mid = 0.5d0 * (z_low + z_high)
f_mid = z_mid * tan(z_mid) - Bix
if (f_mid < 0.0d0) then
z_low = z_mid
else
z_high = z_mid
end if
end do
zeta_x = z_mid
C_x = (4.0d0 * sin(zeta_x)) / (2.0d0 * zeta_x + sin(2.0d0 * zeta_x))
! Solve transcendental root zeta_y * tan(zeta_y) = Biy
z_low = 0.0001d0
z_high = min(1.5707d0, max(0.1d0, atan(Biy)))
if (Biy > 10.0d0) z_high = 1.5700d0
do iter = 1, 60
z_mid = 0.5d0 * (z_low + z_high)
f_mid = z_mid * tan(z_mid) - Biy
if (f_mid < 0.0d0) then
z_low = z_mid
else
z_high = z_mid
end if
end do
zeta_y = z_mid
C_y = (4.0d0 * sin(zeta_y)) / (2.0d0 * zeta_y + sin(2.0d0 * zeta_y))
! Product Solution for 2D Conduction
theta_x0 = C_x * exp(-(zeta_x**2) * Fox)
theta_y0 = C_y * exp(-(zeta_y**2) * Foy)
if (theta_x0 > 1.0d0) theta_x0 = 1.0d0
if (theta_y0 > 1.0d0) theta_y0 = 1.0d0
theta_center = theta_x0 * theta_y0
theta_x_surf = theta_x0 * cos(zeta_x)
theta_y_surf = theta_y0 * cos(zeta_y)
theta_corner = theta_x_surf * theta_y_surf
T_center_C = Tinf_C + (Ti_C - Tinf_C) * theta_center
T_corner_C = Tinf_C + (Ti_C - Tinf_C) * theta_corner
! Total energy removed ratio Q/Qmax
Q_ratio = 1.0d0 - theta_center * (sin(zeta_x)*sin(zeta_y) / (zeta_x*zeta_y))
! Output Results
write(*,'(A)') '============================================================'
write(*,'(A)') ' THERMOFLUIDCALC β 2D TRANSIENT CONDUCTION ENGINE'
write(*,'(A)') '============================================================'
write(*,'(A,F10.1,A,F10.1,A)') 'Billet 2Lx x 2Ly = ', 2.0d0*Lx_mm, ' mm x ', 2.0d0*Ly_mm, ' mm'
write(*,'(A,F10.2,A,E12.4,A)')'Conductivity / Diffusivity= ', k_cond, ' W/m.K / ', alpha_diff, ' m2/s'
write(*,'(A,F10.1,A,F10.1,A)') 'Initial / Ambient Temp = ', Ti_C, ' deg C / ', Tinf_C, ' deg C'
write(*,'(A,F10.1,A,F10.1,A)') 'Quench Time / Convection h= ', time_sec, ' s / ', h_conv, ' W/m2.K'
write(*,'(A)') '------------------------------------------------------------'
write(*,'(A,F10.4,A,F10.4)') 'Biot Numbers (Bix / Biy) = ', Bix, ' / ', Biy
write(*,'(A,F10.4,A,F10.4)') 'Fourier Numbers (Fox/Foy) = ', Fox, ' / ', Foy
write(*,'(A,F10.2,A,F8.4)') 'CENTER TEMPERATURE (T0) = ', T_center_C, ' deg C (theta = ', theta_center, ')'
write(*,'(A,F10.2,A,F8.4)') 'CORNER TEMPERATURE (Tcorn)= ', T_corner_C, ' deg C (theta = ', theta_corner, ')'
write(*,'(A,F10.2,A)') 'Energy Removed (Q/Qmax) = ', Q_ratio * 100.0d0, ' %'
write(*,'(A)') '============================================================'
end program transient_multidim_conduction_adi
π» How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 transient_multidim_conduction_adi.f90 -o transient_multidim_conduction_adi
2. Execution with input.txt redirection:
transient_multidim_conduction_adi < input.txt
π Sample input.txt File Structure
Sample Data:
50.0 80.0 45.0 7850.0 480.0 850.0 30.0 450.0 300.0
Parameter Description:
Half-Width Lx [mm]\nHalf-Height Ly [mm]\nConductivity k [W/(mΒ·K)]\nDensity rho [kg/mΒ³]\nSpecific Heat cp [J/(kgΒ·K)]\nInitial Temp Ti [Β°C]\nAmbient Temp Tinf [Β°C]\nConvection Coeff h [W/(mΒ²Β·K)]\nQuench Time [s]