💻 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.
Stefan Problem — Phase Change
Core Numerical Engine in Fortran 90 • 45 total downloads
stefan_problem.f90
! =========================================================================
! Source File: stefan_problem.f90
! =========================================================================
! ==============================================================================
! Stefan Problem — Phase Change (Solidification/Melting) Calculator
! Neumann Analytical Solution for 1D Semi-Infinite Media
! References: Carslaw & Jaeger Ch.11, Alexiades & Solomon
! ==============================================================================
program stefan_problem
implicit none
! Inputs
integer :: mode ! 1=Solidification, 2=Melting
real(8) :: Tm ! Phase change temp (C)
real(8) :: L ! Latent heat of fusion (kJ/kg)
real(8) :: Ti ! Initial temperature (C)
real(8) :: Ts ! Wall/surface temperature (C)
! Solid properties
real(8) :: ks ! conductivity (W/m-K)
real(8) :: rhos ! density (kg/m3)
real(8) :: Cps ! specific heat (J/kg-K)
! Liquid properties
real(8) :: kl ! conductivity (W/m-K)
real(8) :: rhol ! density (kg/m3)
real(8) :: Cpl ! specific heat (J/kg-K)
! Dimensions & simulation time
real(8) :: d_target ! target thickness for phase change (mm)
real(8) :: x_eval ! evaluation depth (mm)
real(8) :: t_sim ! elapsed time (s)
! Derived properties
real(8) :: alphas ! solid diffusivity (m2/s)
real(8) :: alphal ! liquid diffusivity (m2/s)
real(8) :: Ste ! Stefan number
real(8) :: lambda ! similarity growth constant
! Outputs
real(8) :: front_pos ! front position s(t) (m)
real(8) :: time_to_d ! solidification/melting time for d (s)
real(8) :: temp_at_x ! temperature at x_eval (C)
real(8) :: energy_transferred ! energy released/absorbed (kJ/m2)
! Constants & loop variables
real(8) :: pi, eps
integer :: iter
real(8) :: a, b, c, fa, fc
character(len=20) :: mode_str
pi = 4.0d0 * datan(1.0d0)
eps = 1.0d-12
! Read standard inputs
read(*, *, iostat=iter) mode
if (iter /= 0) then
print *, "Error: Invalid input format (mode)."
stop
endif
read(*, *) Tm
read(*, *) L
read(*, *) Ti
read(*, *) Ts
read(*, *) ks, rhos, Cps
read(*, *) kl, rhol, Cpl
read(*, *) d_target
read(*, *) x_eval
read(*, *) t_sim
! Sanity Checks & Input Validation
if (L <= 0.0d0 .or. ks <= 0.0d0 .or. rhos <= 0.0d0 .or. Cps <= 0.0d0 &
.or. kl <= 0.0d0 .or. rhol <= 0.0d0 .or. Cpl <= 0.0d0) then
print *, "Error: Physical properties must be positive."
stop
endif
if (d_target < 0.0d0 .or. x_eval < 0.0d0 .or. t_sim <= 0.0d0) then
print *, "Error: Target thickness, depth, and time must be non-negative."
stop
endif
if (mode == 1) then
mode_str = "Solidification"
if (Ts >= Tm) then
print *, "Error: For Solidification, wall temperature (Ts) must be below melting temperature (Tm)."
stop
endif
if (Ti < Tm) then
print *, "Error: For Solidification, initial temperature (Ti) must be above or at melting temperature (Tm)."
stop
endif
else if (mode == 2) then
mode_str = "Melting"
if (Ts <= Tm) then
print *, "Error: For Melting, wall temperature (Ts) must be above melting temperature (Tm)."
stop
endif
if (Ti > Tm) then
print *, "Error: For Melting, initial temperature (Ti) must be below or at melting temperature (Tm)."
stop
endif
else
print *, "Error: Invalid Mode (must be 1 for Solidification or 2 for Melting)."
stop
endif
! Diffusivities
alphas = ks / (rhos * Cps)
alphal = kl / (rhol * Cpl)
! Stefan Number
if (mode == 1) then
Ste = Cps * (Tm - Ts) / (L * 1000.0d0)
else
Ste = Cpl * (Ts - Tm) / (L * 1000.0d0)
endif
! Bisection Solver for Lambda
a = 1.0d-8
b = 10.0d0
fa = evaluate_f(a)
if (fa * evaluate_f(b) > 0.0d0) then
print *, "Error: Root is not bracketed in [10^-8, 10.0]."
stop
endif
do iter = 1, 100
c = (a + b) / 2.0d0
fc = evaluate_f(c)
if (abs(fc) < eps .or. (b - a)/2.0d0 < eps) then
lambda = c
exit
endif
if (fa * fc < 0.0d0) then
b = c
else
a = c
fa = fc
endif
enddo
! Calculations
if (mode == 1) then
! Growing phase is solid
front_pos = 2.0d0 * lambda * sqrt(alphas * t_sim)
time_to_d = (d_target / 1000.0d0)**2 / (4.0d0 * lambda**2 * alphas)
! Temperature at x_eval
if (x_eval / 1000.0d0 < front_pos) then
! Solid region
temp_at_x = Ts + (Tm - Ts) * erf( (x_eval / 1000.0d0) / (2.0d0 * sqrt(alphas * t_sim)) ) / erf(lambda)
else
! Liquid region
temp_at_x = Ti - (Ti - Tm) * erfc( (x_eval / 1000.0d0) / (2.0d0 * sqrt(alphal * t_sim)) ) / &
erfc(lambda * sqrt(alphas / alphal))
endif
! Energy released (kJ/m2)
energy_transferred = 2.0d0 * ks * (Tm - Ts) / erf(lambda) * sqrt(t_sim / (pi * alphas)) / 1000.0d0
else
! Growing phase is liquid
front_pos = 2.0d0 * lambda * sqrt(alphal * t_sim)
time_to_d = (d_target / 1000.0d0)**2 / (4.0d0 * lambda**2 * alphal)
! Temperature at x_eval
if (x_eval / 1000.0d0 < front_pos) then
! Liquid region
temp_at_x = Ts - (Ts - Tm) * erf( (x_eval / 1000.0d0) / (2.0d0 * sqrt(alphal * t_sim)) ) / erf(lambda)
else
! Solid region
temp_at_x = Ti + (Tm - Ti) * erfc( (x_eval / 1000.0d0) / (2.0d0 * sqrt(alphas * t_sim)) ) / &
erfc(lambda * sqrt(alphal / alphas))
endif
! Energy absorbed (kJ/m2)
energy_transferred = 2.0d0 * kl * (Ts - Tm) / erf(lambda) * sqrt(t_sim / (pi * alphal)) / 1000.0d0
endif
! ==========================================================================
! Print Results in a beautiful engineering report format
! ==========================================================================
print *, "=========================================================================="
print *, " STEFAN PROBLEM PHASE CHANGE REPORT (NEUMANN SOLUTION) "
print *, "=========================================================================="
print "(A, A20)", " Operation Mode: ", trim(mode_str)
print "(A, F10.2, A)", " Phase Change Temp (Tm): ", Tm, " C"
print "(A, F10.2, A)", " Latent Heat of Fusion (L):", L, " kJ/kg"
print "(A, F10.2, A)", " Initial Temp (Ti): ", Ti, " C"
print "(A, F10.2, A)", " Wall/Surface Temp (Ts): ", Ts, " C"
print *
print *, "--------------------------------------------------------------------------"
print *, " MATERIAL PROPERTIES "
print *, "--------------------------------------------------------------------------"
print *, " Parameter Solid Phase Liquid Phase"
print "(A, F18.3, F21.3)", " Conductivity (k): ", ks, kl
print "(A, F18.1, F21.1)", " Density (rho): ", rhos, rhol
print "(A, F18.1, F21.1)", " Spec. Heat (Cp): ", Cps, Cpl
print "(A, ES18.4, ES21.4)", " Diffusivity (alpha):", alphas, alphal
print *
print *, "--------------------------------------------------------------------------"
print *, " COMPUTATION & DIMENSIONLESS METRICS "
print *, "--------------------------------------------------------------------------"
print "(A, F12.5)", " Stefan Number (Ste): ", Ste
print "(A, F12.6)", " Growth Eigenvalue (lambda): ", lambda
print *
print *, "--------------------------------------------------------------------------"
print *, " DETAILED PERFORMANCE RESULTS "
print *, "--------------------------------------------------------------------------"
print "(A, F12.2, A)", " Elapsed Simulation Time: ", t_sim, " s"
print "(A, F12.2, A)", " Phase Interface Position s(t): ", front_pos * 1000.0d0, " mm"
print "(A, F8.2, A, F12.2, A)", " Phase Change Time for d =", d_target, " mm: ", time_to_d, " s"
print "(A, F8.2, A, F12.2, A)", " Temperature at depth x =", x_eval, " mm: ", temp_at_x, " C"
if (mode == 1) then
print "(A, F12.2, A)", " Total Thermal Energy Released: ", energy_transferred, " kJ/m2"
else
print "(A, F12.2, A)", " Total Thermal Energy Absorbed: ", energy_transferred, " kJ/m2"
endif
print *, "=========================================================================="
print *, " References:"
print *, " 1. Carslaw, H. S. & Jaeger, J. C., Conduction of Heat in Solids, Ch. 11"
print *, " 2. Alexiades, V. & Solomon, A. D., Mathematical Modeling of Phase Change"
print *, "=========================================================================="
contains
! Helper function to evaluate the transcendental equation residual f(lambda)
real(8) function evaluate_f(l)
real(8), intent(in) :: l
real(8) :: term1, term2, y
! Term 1: Ste / (exp(l^2) * erf(l))
if (l > 20.0d0) then
term1 = 0.0d0
else
term1 = Ste / (exp(l**2) * erf(l))
endif
! Term 2: solid/liquid coupling term
if (mode == 1) then
! Solidification
y = l * sqrt(alphas / alphal)
term2 = (kl * sqrt(alphas) * (Ti - Tm)) / &
(ks * sqrt(alphal) * (Tm - Ts)) * Ste * exp_erfc_reciprocal(y)
else
! Melting
y = l * sqrt(alphal / alphas)
term2 = (ks * sqrt(alphal) * (Tm - Ti)) / &
(kl * sqrt(alphas) * (Ts - Tm)) * Ste * exp_erfc_reciprocal(y)
endif
evaluate_f = term1 - term2 - l * sqrt(pi)
end function evaluate_f
! Helper to safely evaluate 1 / (exp(y^2) * erfc(y)) without overflow
real(8) function exp_erfc_reciprocal(y)
real(8), intent(in) :: y
if (y < 10.0d0) then
exp_erfc_reciprocal = 1.0d0 / (exp(y**2) * erfc(y))
else
! Asymptotic expansion for y -> infinity
exp_erfc_reciprocal = (y * sqrt(pi)) / &
(1.0d0 - 0.5d0 / (y**2) + 0.75d0 / (y**4))
endif
end function exp_erfc_reciprocal
end program stefan_problem
Solver Description
Calculate moving boundaries, temperature distributions, and heat rates in 1D phase change solidification or melting using the Neumann analytical solution.
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 stefan_problem.f90 -o stefan_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
stefan_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Stefan Mode (1=Freezing/Solidification, 2=Melting/Liquefaction)
1
! Melting Temperature Tm [°C]
0.0
! Latent Heat of Fusion L [kJ/kg]
333.5
! Initial Liquid Temp Ti [°C]
10.0
! Applied Surface Temp Ts [°C]
-15.0
! Solid phase thermal conductivity ks [W/m-K]
2.22
! Solid phase density rhos [kg/m3]
917.0
! Solid phase specific heat Cps [J/kg-K]
2050.0
! Liquid phase thermal conductivity kl [W/m-K]
0.57
! Liquid phase density rhol [kg/m3]
1000.0
! Liquid phase specific heat Cpl [J/kg-K]
4184.0
! Target freeze depth [mm]
50.0
! Evaluation depth x_eval [mm]
10.0
! Simulation time t [s]
3600.0
1
! Melting Temperature Tm [°C]
0.0
! Latent Heat of Fusion L [kJ/kg]
333.5
! Initial Liquid Temp Ti [°C]
10.0
! Applied Surface Temp Ts [°C]
-15.0
! Solid phase thermal conductivity ks [W/m-K]
2.22
! Solid phase density rhos [kg/m3]
917.0
! Solid phase specific heat Cps [J/kg-K]
2050.0
! Liquid phase thermal conductivity kl [W/m-K]
0.57
! Liquid phase density rhol [kg/m3]
1000.0
! Liquid phase specific heat Cpl [J/kg-K]
4184.0
! Target freeze depth [mm]
50.0
! Evaluation depth x_eval [mm]
10.0
! Simulation time t [s]
3600.0