💻 Fortran Source Code Library
Run calculations locally on your own machine. View code structure, read technical explanations, and download compilation packages including sample input files.
Boiling & Condensation Phase-Change Solver
Core Numerical Engine in Fortran 90 • 0 total downloads
boiling_condensation.f90
! =========================================================================
! Source File: boiling_condensation.f90
! =========================================================================
program boiling_condensation
implicit none
! Inputs
integer :: mode ! 1 = Pool Boiling, 2 = Film Condensation
integer :: geom_type ! For Cond: 1=Vertical Plate, 2=Cylinder, 3=Sphere. For Boiling: ignored.
double precision :: dim1, dim2
double precision :: T_s, T_sat
double precision :: rho_l, mu_l, k_l, Pr_l, Cp_l
double precision :: rho_v, mu_v, k_v, Pr_v, Cp_v
double precision :: h_fg, sigma, C_sf, emiss
! Constants
double precision, parameter :: g = 9.80665d0
double precision, parameter :: sigma_sb = 5.67d-8
double precision, parameter :: PI = 3.141592653589793d0
! Solved variables
double precision :: dT_e, dT, h_fg_mod, h_avg, Q, A_s, flux, q_max, q_min, L_c
double precision :: h_film, h_rad, T_s_K, T_sat_K
double precision :: mass_flow
! Local profile
integer :: i, n_points
double precision :: x, dx, local_dT_e, local_flux, log_dT, local_delta
character(len=50) :: regime_name
integer :: iostat_val
! Read inputs
read(*,*,iostat=iostat_val) mode
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid mode input.'
stop
end if
read(*,*,iostat=iostat_val) geom_type
read(*,*,iostat=iostat_val) dim1
read(*,*,iostat=iostat_val) dim2
read(*,*,iostat=iostat_val) T_s
read(*,*,iostat=iostat_val) T_sat
read(*,*,iostat=iostat_val) rho_l
read(*,*,iostat=iostat_val) mu_l
read(*,*,iostat=iostat_val) k_l
read(*,*,iostat=iostat_val) Pr_l
read(*,*,iostat=iostat_val) Cp_l
read(*,*,iostat=iostat_val) rho_v
read(*,*,iostat=iostat_val) mu_v
read(*,*,iostat=iostat_val) k_v
read(*,*,iostat=iostat_val) Pr_v
read(*,*,iostat=iostat_val) Cp_v
read(*,*,iostat=iostat_val) h_fg
read(*,*,iostat=iostat_val) sigma
read(*,*,iostat=iostat_val) C_sf
read(*,*,iostat=iostat_val) emiss
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all simulation parameters.'
stop
end if
! Basic validations
if (dim1 <= 0.0d0) dim1 = 0.01d0 ! Safe cylinder diameter default
if (rho_l <= 0.0d0 .or. h_fg <= 0.0d0 .or. sigma <= 0.0d0) then
write(*,*) 'ERROR: Density, latent heat, and surface tension must be positive.'
stop
end if
if (mode == 1) then
! ============================================================
! POOL BOILING ANALYSIS
! ============================================================
dT_e = T_s - T_sat
if (dT_e < 0.0d0) then
dT_e = 0.0d0
end if
! Zuber CHF calculation
q_max = 0.149d0 * h_fg * rho_v * ((g * sigma * (rho_l - rho_v)) / (rho_v**2))**0.25d0
! Leidenfrost point minimum heat flux
q_min = 0.09d0 * rho_v * h_fg * ((g * sigma * (rho_l - rho_v)) / ((rho_l + rho_v)**2))**0.25d0
! Saturated pool boiling curve piecewise modeling
T_s_K = T_s + 273.15d0
T_sat_K = T_sat + 273.15d0
! Nucleate boiling Rohsenow correlation
if (dT_e > 0.0d0) then
flux = mu_l * h_fg * sqrt((g * (rho_l - rho_v)) / sigma) * &
((Cp_l * dT_e) / (C_sf * h_fg * Pr_l))**3
else
flux = 0.0d0
end if
! Check Bromley film boiling (at Leidenfrost or higher)
h_fg_mod = h_fg + 0.8d0 * Cp_v * dT_e
h_film = 0.62d0 * ((g * rho_v * (rho_l - rho_v) * k_v**3 * h_fg_mod) / &
(mu_v * dim1 * max(dT_e, 1.0d-5)))**0.25d0
h_rad = (emiss * sigma_sb * (T_s_K**4 - T_sat_K**4)) / max(dT_e, 1.0d-5)
h_avg = h_film + 0.75d0 * h_rad
! Piecewise heat flux selection matching Nukiyama regimes
if (dT_e < 5.0d0) then
regime_name = "Natural Convection Boiling"
flux = 500.0d0 * dT_e**1.25d0
else if (dT_e < 30.0d0) then
regime_name = "Nucleate Boiling"
! Rohsenow flux matches standard nucleate regime
if (flux > q_max) flux = q_max
else if (dT_e < 120.0d0) then
regime_name = "Transition Boiling"
! Log-log interpolation between CHF and Leidenfrost points
log_dT = (log10(dT_e) - log10(30.0d0)) / (log10(120.0d0) - log10(30.0d0))
flux = 10.0d0**(log10(q_max) + log_dT * (log10(q_min) - log10(q_max)))
else
regime_name = "Film Boiling"
flux = h_avg * dT_e
if (flux < q_min) flux = q_min
end if
! Output Boiling results
write(*,'(A)') '============================================================'
write(*,'(A)') ' POOL BOILING HEAT TRANSFER CURVE ANALYSIS'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A,A)') ' Boiling Regime = ', trim(regime_name)
write(*,'(A,F12.2,A)') ' Excess Temperature (dTe)= ', dT_e, ' deg-C'
write(*,'(A,F12.2,A)') ' Surface Temperature (Ts)= ', T_s, ' deg-C'
write(*,'(A,F12.2,A)') ' Saturation Temp (Tsat) = ', T_sat, ' deg-C'
write(*,*)
write(*,'(A)') '--- TRANSITION BOUNDARIES -----------------------------------'
write(*,'(A,ES12.4,A)') ' Critical Heat Flux (CHF)= ', q_max, ' W/m2'
write(*,'(A,ES12.4,A)') ' Leidenfrost Flux (min) = ', q_min, ' W/m2'
write(*,*)
write(*,'(A)') '--- HEAT TRANSFER RESULTS -----------------------------------'
write(*,'(A,ES12.4,A)') ' Calculated Heat Flux (q")= ', flux, ' W/m2'
write(*,'(A,F12.4,A)') ' Film HTC (Convective) = ', h_film, ' W/m2.K'
write(*,'(A,F12.4,A)') ' Film HTC (Radiation) = ', h_rad, ' W/m2.K'
write(*,'(A,F12.4,A)') ' Total HTC (h_total) = ', h_avg, ' W/m2.K'
write(*,*)
! Generate log boiling curve dataset
write(*,'(A)') '--- LOCAL PROFILE ALONG SURFACE ----------------------------'
write(*,'(A)') ' dTe [C] Heat_Flux [W/m2]'
write(*,'(A)') ' ----------------------------------------------------------'
n_points = 40
do i = 0, n_points
log_dT = 3.0d0 * dble(i) / dble(n_points) ! dT_e from 10^0 (1) to 10^3 (1000)
local_dT_e = 10.0d0**log_dT
! Piecewise calculation for profile
if (local_dT_e < 5.0d0) then
local_flux = 500.0d0 * local_dT_e**1.25d0
else if (local_dT_e < 30.0d0) then
local_flux = mu_l * h_fg * sqrt((g * (rho_l - rho_v)) / sigma) * &
((Cp_l * local_dT_e) / (C_sf * h_fg * Pr_l))**3
if (local_flux > q_max) local_flux = q_max
else if (local_dT_e < 120.0d0) then
log_dT = (log10(local_dT_e) - log10(30.0d0)) / (log10(120.0d0) - log10(30.0d0))
local_flux = 10.0d0**(log10(q_max) + log_dT * (log10(q_min) - log10(q_max)))
else
h_fg_mod = h_fg + 0.8d0 * Cp_v * local_dT_e
h_film = 0.62d0 * ((g * rho_v * (rho_l - rho_v) * k_v**3 * h_fg_mod) / &
(mu_v * dim1 * local_dT_e))**0.25d0
h_rad = (emiss * sigma_sb * (((T_sat + local_dT_e + 273.15d0)**4) - &
(T_sat + 273.15d0)**4)) / local_dT_e
local_flux = (h_film + 0.75d0 * h_rad) * local_dT_e
if (local_flux < q_min) local_flux = q_min
end if
write(*,'(F10.4,4X,ES14.6)') local_dT_e, local_flux
end do
else
! ============================================================
! FILM CONDENSATION ANALYSIS
! ============================================================
dT = T_sat - T_s
if (dT <= 0.0d0) then
write(*,*) 'ERROR: Surface must be colder than saturation temp for condensation.'
stop
end if
h_fg_mod = h_fg + 0.68d0 * Cp_l * dT
select case (geom_type)
case (1)
regime_name = "Vertical Plate Condensation"
L_c = dim1
if (dim2 <= 0.0d0) dim2 = 1.0d0 ! width default
A_s = dim1 * dim2
h_avg = 0.943d0 * ((g * rho_l * (rho_l - rho_v) * k_l**3 * h_fg_mod) / &
(mu_l * L_c * dT))**0.25d0
case (2)
regime_name = "Horizontal Cylinder Condensation"
L_c = dim1 ! Diameter
if (dim2 <= 0.0d0) dim2 = 1.0d0 ! cylinder length
A_s = PI * dim1 * dim2
h_avg = 0.729d0 * ((g * rho_l * (rho_l - rho_v) * k_l**3 * h_fg_mod) / &
(mu_l * L_c * dT))**0.25d0
case (3)
regime_name = "Sphere Condensation"
L_c = dim1
A_s = PI * dim1**2
h_avg = 0.815d0 * ((g * rho_l * (rho_l - rho_v) * k_l**3 * h_fg_mod) / &
(mu_l * L_c * dT))**0.25d0
case default
write(*,*) 'ERROR: Invalid geometry type for condensation.'
stop
end select
Q = h_avg * A_s * dT
mass_flow = Q / h_fg_mod
! Output Condensation results
write(*,'(A)') '============================================================'
write(*,'(A)') ' FILM CONDENSATION HEAT TRANSFER SIZING'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A,A)') ' Condensation Regime = ', trim(regime_name)
write(*,'(A,F12.4,A)') ' Characteristic Length = ', L_c, ' m'
write(*,'(A,F12.4,A)') ' Heat Transfer Area (As) = ', A_s, ' m2'
write(*,'(A,F12.2,A)') ' Saturation Temp (Tsat) = ', T_sat, ' deg-C'
write(*,'(A,F12.2,A)') ' Surface Temperature (Ts)= ', T_s, ' deg-C'
write(*,'(A,F12.2,A)') ' Temp Difference (dT) = ', dT, ' deg-C'
write(*,*)
write(*,'(A)') '--- HEAT TRANSFER RESULTS -----------------------------------'
write(*,'(A,F12.4,A)') ' Modified Latent Heat = ', h_fg_mod / 1000.0d0, ' kJ/kg'
write(*,'(A,F12.4,A)') ' Average Condensation h = ', h_avg, ' W/m2.K'
write(*,'(A,F12.4,A)') ' Total Heat Transfer (Q) = ', Q, ' W'
write(*,'(A,ES12.4,A)') ' Condensation Mass Rate = ', mass_flow, ' kg/s'
write(*,*)
! Generate local boundary layer profile data (film thickness delta along x)
write(*,'(A)') '--- LOCAL PROFILE ALONG SURFACE ----------------------------'
write(*,'(A)') ' x [m] film_thickness [mm]'
write(*,'(A)') ' ----------------------------------------------------------'
n_points = 40
dx = L_c / dble(n_points)
do i = 1, n_points
x = dble(i) * dx
! Nusselt film thickness relation
local_delta = ((4.0d0 * mu_l * k_l * dT * x) / (g * rho_l * (rho_l - rho_v) * h_fg_mod))**0.25d0
write(*,'(F8.4,4X,F12.6)') x, local_delta * 1000.0d0
end do
end if
end program boiling_condensation
Solver Description
Solves phase-change heat transfer. For pool boiling, it determines the regime along the Nukiyama curve (nucleate, critical heat flux using Zuber's limit, or film boiling using Rohsenow and Bromley equations). For film condensation, it uses Nusselt's analytical liquid film drainage theory.
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 -ffree-line-length-none boiling_condensation.f90 -o boiling_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
boiling_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Mode (1=Pool Boiling, 2=Condensation)
1
! Surface Temp ($T_s$) [°C]
115.0
! Saturation Temp ($T_{sat}$) [°C]
100.0
! Liquid Density ($\rho_l$) [kg/m³]
957.0
! Liquid Viscosity ($\mu_l$) [Pa-s]
0.00028
! Liquid Conductivity ($k_l$) [W/m-K]
0.68
! Liquid Specific Heat ($C_{p,l}$) [J/kg-K]
4217.0
! Latent Heat of Vaporization ($h_{fg}$) [J/kg]
2257000.0
! Surface Tension ($\sigma$) [N/m] (ignored for Condensation)
0.0589
! Surface Coefficient ($C_{sf}$) [0-1] (ignored for Condensation)
0.013
! Characteristic Dimension (Diameter or Length) [m]
0.05
1
! Surface Temp ($T_s$) [°C]
115.0
! Saturation Temp ($T_{sat}$) [°C]
100.0
! Liquid Density ($\rho_l$) [kg/m³]
957.0
! Liquid Viscosity ($\mu_l$) [Pa-s]
0.00028
! Liquid Conductivity ($k_l$) [W/m-K]
0.68
! Liquid Specific Heat ($C_{p,l}$) [J/kg-K]
4217.0
! Latent Heat of Vaporization ($h_{fg}$) [J/kg]
2257000.0
! Surface Tension ($\sigma$) [N/m] (ignored for Condensation)
0.0589
! Surface Coefficient ($C_{sf}$) [0-1] (ignored for Condensation)
0.013
! Characteristic Dimension (Diameter or Length) [m]
0.05