Stefan Moving Boundary Solidification Front

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

📦 Download Complete ZIP 📥 Download .f90 Code 📄 Download input.txt

🔬

Solver Purpose & Physical Scope

Solve two-phase moving boundary Stefan phase change problems: solid front position s(t), total freezing time, solid Stefan number (Stes), and instantaneous cooling heat flux.

📂 Discipline: Conduction ⚡ Precision: IEEE-754 64-bit Real(`real(8)`) 📥 Total Downloads: 401 times 📄 Source File: stefan_phase_change_moving_boundary.f90

📁 calcul/Conduction / stefan_phase_change_moving_boundary.f90

program stefan_phase_change_moving_boundary
    implicit none
    integer :: iostat_val, iter
    double precision :: L_target_mm, Tm_melt_C, T0_wall_C, Ti_liquid_C, Lf_latent_kJkg
    double precision :: ks_solid, rhos_solid, cps_solid, kl_liquid, rhol_liquid, cpl_liquid
    double precision :: L_m, alpha_s, alpha_l, Ste_s, Ste_l, lambda_val
    double precision :: lam_low, lam_high, lam_mid, f_mid, t_freeze_sec, t_freeze_hr
    double precision :: s_front_1hr_mm, q_flux_init_kWm2
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) L_target_mm     ! Target Thickness to Freeze [mm] (e.g. 50.0)
    read(*,*,iostat=iostat_val) Tm_melt_C       ! Melting / Freezing Temp [deg C] (e.g. 0.0 for water, 58.0 for paraffin)
    read(*,*,iostat=iostat_val) T0_wall_C       ! Cold Wall Temp [deg C] (e.g. -15.0)
    read(*,*,iostat=iostat_val) Ti_liquid_C     ! Initial Liquid Temp [deg C] (e.g. 10.0)
    read(*,*,iostat=iostat_val) Lf_latent_kJkg  ! Latent Heat of Fusion [kJ/kg] (e.g. 333.5 for ice)
    read(*,*,iostat=iostat_val) ks_solid        ! Solid Conductivity [W/(m.K)] (e.g. 2.22 for ice)
    read(*,*,iostat=iostat_val) rhos_solid      ! Solid Density [kg/m3] (e.g. 917.0)
    read(*,*,iostat=iostat_val) cps_solid       ! Solid Specific Heat [J/(kg.K)] (e.g. 2050.0)
    read(*,*,iostat=iostat_val) kl_liquid       ! Liquid Conductivity [W/(m.K)] (e.g. 0.58)
    read(*,*,iostat=iostat_val) rhol_liquid     ! Liquid Density [kg/m3] (e.g. 1000.0)
    read(*,*,iostat=iostat_val) cpl_liquid      ! Liquid Specific Heat [J/(kg.K)] (e.g. 4180.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for Stefan phase change calculation.'
        stop
    end if

    if (L_target_mm <= 0.0d0 .or. Lf_latent_kJkg <= 0.0d0 .or. ks_solid <= 0.0d0) then
        write(*,*) 'ERROR: Thickness, latent heat, and conductivity must be positive.'
        stop
    end if

    L_m = L_target_mm * 1.0d-3
    alpha_s = ks_solid / (rhos_solid * cps_solid)
    alpha_l = kl_liquid / (rhol_liquid * cpl_liquid)

    Ste_s = (cps_solid * (Tm_melt_C - T0_wall_C)) / (Lf_latent_kJkg * 1000.0d0)
    Ste_l = (cpl_liquid * (Ti_liquid_C - Tm_melt_C)) / (Lf_latent_kJkg * 1000.0d0)

    ! Solve transcendental root for lambda: lambda * exp(lambda^2) * erf(lambda) = Ste_s / sqrt(PI)
    lam_low = 0.001d0
    lam_high = 3.0d0
    do iter = 1, 60
        lam_mid = 0.5d0 * (lam_low + lam_high)
        ! Approximate erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10) for small x or erf intrinsic
        f_mid = lam_mid * exp(lam_mid**2) * erf(lam_mid) - (Ste_s / sqrt(PI))
        if (f_mid < 0.0d0) then
            lam_low = lam_mid
        else
            lam_high = lam_mid
        end if
    end do
    lambda_val = lam_mid

    ! Freezing Time for Target Thickness
    t_freeze_sec = (L_m**2) / (4.0d0 * (lambda_val**2) * alpha_s)
    t_freeze_hr = t_freeze_sec / 3600.0d0

    ! Front position after 1 hour [mm]
    s_front_1hr_mm = 2.0d0 * lambda_val * sqrt(alpha_s * 3600.0d0) * 1000.0d0

    ! Heat flux at t = 60s
    q_flux_init_kWm2 = (ks_solid * (Tm_melt_C - T0_wall_C) / (sqrt(PI * alpha_s * 60.0d0) * erf(lambda_val))) / 1000.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — TWO-PHASE STEFAN MOVING BOUNDARY ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Target Freeze Thickness   = ', L_target_mm, ' mm (Wall T0 = ', T0_wall_C, ' C)'
    write(*,'(A,F10.1,A,F10.1,A)') 'Melting Tm / Initial Liq  = ', Tm_melt_C, ' C / ', Ti_liquid_C, ' C'
    write(*,'(A,F10.4,A,F10.4)')   'Stefan Number Ste_s / Ste_l= ', Ste_s, ' / ', Ste_l
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.4)')           'FRONT GROWTH CONSTANT lam = ', lambda_val
    write(*,'(A,F10.2,A,F10.2,A)') 'TOTAL FREEZING TIME       = ', t_freeze_hr, ' hours (', t_freeze_sec, ' s)'
    write(*,'(A,F10.2,A)')  'Front Thickness after 1 hr= ', s_front_1hr_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Instant Heat Flux (t=60s) = ', q_flux_init_kWm2, ' kW/m2'
    write(*,'(A)') '============================================================'

end program stefan_phase_change_moving_boundary

💻 How to Compile & Run Locally

1. Compilation (GNU Fortran / Intel oneAPI):

gfortran -O3 stefan_phase_change_moving_boundary.f90 -o stefan_phase_change_moving_boundary

2. Execution with input.txt redirection:

stefan_phase_change_moving_boundary < input.txt

📄 Sample input.txt File Structure

Sample Data:

Parameter Description:

Target Thickness L [mm]\nMelting Temp Tm [°C]\nCold Wall Temp T0 [°C]\nInitial Liquid Temp Ti [°C]\nLatent Heat Lf [kJ/kg]\nSolid Conductivity ks [W/(m·K)]\nSolid Density rhos [kg/m³]\nSolid Specific Heat cps [J/(kg·K)]\nLiquid Conductivity kl [W/(m·K)]\nLiquid Density rhol [kg/m³]\nLiquid Specific Heat cpl [J/(kg·K)]