🔬
Solver Purpose & Physical Scope
Phase-change heat transfer includes boiling and condensation. Nucleate pool boiling is modeled using Rohsenow's correlation for heat flux:
📂 Discipline: Convection
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 277 times
📄 Source File:
boiling_condensation.f90
📁 calcul/Convection /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 boiling_condensation.f90 -o boiling_calc
2. Execution with input.txt redirection:
boiling_calc < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 1 0.025 1.0 110.0 100.0 958.0 2.82e-4 0.68 1.75 4217.0 0.598 1.2e-5 0.025 1.0 2030.0 2257000.0 0.0589 0.013 0.8
Parameter Description:
Mode (1=Pool Boiling, 2=Film Condensation) Geometry (1=Horizontal Wire/Tube, 2=Flat Plate, 3=Sphere) Diameter/Height [m] Length/Width [m] Surface Temp Ts [°C] Saturation Temp Tsat [°C] Liquid Density rhol [kg/m3] Liquid Viscosity mul [Pa-s] Liquid Conductivity kl [W/m-K] Liquid Prandtl Prl Liquid Specific Heat Cpl [J/kg-K] Vapor Density rhov [kg/m3] Vapor Viscosity muv [Pa-s] Vapor Conductivity kv [W/m-K] Vapor Prandtl Prv Vapor Specific Heat Cpv [J/kg-K] Latent Heat hfg [J/kg] Surface Tension sigma [N/m] Boiling correlation constant Csf (dummy) Surface Emissivity