⚑ Fortran 90 / 2008 Double Precision
πŸ“₯ 321 Downloads

2D Transient Conduction & Heisler Multi-D

Standalone, self-contained numerical routine. Verify algorithms, inspect boundary condition equations, or compile locally for batch parametric runs.

πŸ”¬

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]