🔬
Solver Purpose & Physical Scope
Calculates droplet Sauter Mean Diameter (SMD, D32) using Lefebvre-Hiroyasu hybrid models, identifies liquid breakup regimes (Rayleigh, Wind-Induced, Catastrophic Atomization) via Weber and Ohnesorge dimensionless numbers, computes spray cone angles and intact liquid core length, and generates Rosin-Rammler droplet volume distributions.
📂 Discipline: Fluid
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 211 times
📄 Source File:
spray_atomization.f90
📁 calcul/FluidMechanics /
spray_atomization.f90
program spray_atomization
implicit none
integer :: i, iostat_val, n_points
double precision :: d_nozzle_mm, d_m, deltaP_bar, deltaP_Pa, Cd
double precision :: rho_L, mu_L, sigma, rho_g
double precision :: A_orifice, V_ideal, V_inj, mdot_L, Q_L
double precision :: Re_L, We_L, We_g, Oh
double precision :: term1, term2, smd_m, D32, D10, Dv50
double precision :: rho_ratio, theta_rad, spray_angle, Lb_m, breakup_length
double precision :: q_rr, X_rr, d_val, cum_vol, pdf_val
double precision :: p_var, p_var_Pa, m_var, smd_var
double precision, parameter :: PI = 3.141592653589793d0
character(len=64) :: regime_name
! Read inputs from standard input
read(*,*,iostat=iostat_val) d_nozzle_mm
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid nozzle diameter input.'
stop
end if
read(*,*,iostat=iostat_val) deltaP_bar
read(*,*,iostat=iostat_val) Cd
read(*,*,iostat=iostat_val) rho_L
read(*,*,iostat=iostat_val) mu_L
read(*,*,iostat=iostat_val) sigma
read(*,*,iostat=iostat_val) rho_g
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all spray atomization inputs.'
stop
end if
! Validation
if (d_nozzle_mm <= 0.0d0 .or. deltaP_bar <= 0.0d0 .or. Cd <= 0.0d0) then
write(*,*) 'ERROR: Geometric and pressure inputs must be positive.'
stop
end if
if (rho_L <= 0.0d0 .or. mu_L <= 0.0d0 .or. sigma <= 0.0d0 .or. rho_g <= 0.0d0) then
write(*,*) 'ERROR: Physical fluid properties must be positive.'
stop
end if
! Conversions
d_m = d_nozzle_mm * 1.0d-3
deltaP_Pa = deltaP_bar * 1.0d5
A_orifice = (PI / 4.0d0) * (d_m**2)
! Hydrodynamics & Flow Rates
V_ideal = sqrt(2.0d0 * deltaP_Pa / rho_L)
V_inj = Cd * V_ideal
mdot_L = Cd * A_orifice * sqrt(2.0d0 * rho_L * deltaP_Pa)
Q_L = mdot_L / rho_L
! Non-dimensional numbers
Re_L = (rho_L * V_inj * d_m) / mu_L
We_L = (rho_L * (V_inj**2) * d_m) / sigma
We_g = (rho_g * (V_inj**2) * d_m) / sigma
Oh = mu_L / sqrt(rho_L * sigma * d_m)
! Breakup Regime Identification
if (Oh < 0.1d0 .and. We_L < 10.0d0) then
regime_name = 'Rayleigh Capillary Breakup'
else if (We_g < 13.0d0) then
regime_name = 'First Wind-Induced Regime'
else if (We_g < 40.3d0) then
regime_name = 'Second Wind-Induced Regime'
else
regime_name = 'Catastrophic Shear Atomization'
end if
! Sauter Mean Diameter (SMD D32) — Lefebvre-Hiroyasu Hybrid Model
term1 = 2.25d0 * (sigma**0.25d0) * (mu_L**0.16d0) * (rho_L**0.20d0)
term2 = (mdot_L**0.22d0) * (deltaP_Pa**(-0.43d0))
smd_m = term1 * term2
if (We_g > 40.0d0) then
smd_m = smd_m * ((40.0d0 / We_g)**0.12d0)
end if
D32 = max(1.0d0, smd_m * 1.0d6) ! in micrometers
D10 = D32 * 0.65d0
Dv50 = D32 * 1.18d0 ! Mass Median Diameter
! Spray Cone Angle (Arai/Hiroyasu)
rho_ratio = rho_g / rho_L
theta_rad = 0.057d0 * (rho_ratio**0.10d0) * (max(100.0d0, Re_L)**0.25d0)
spray_angle = min(120.0d0, max(5.0d0, theta_rad * (180.0d0 / PI) * 1.8d0))
! Liquid Core Length
Lb_m = 7.0d0 * d_m * sqrt(rho_L / rho_g) * (1.0d0 + 0.4d0 * Oh)
breakup_length = Lb_m * 1.0d3 ! in mm
! Rosin-Rammler parameters
q_rr = 2.85d0
X_rr = Dv50 / (log(2.0d0)**(1.0d0 / q_rr))
! Output Results
write(*,'(A)') '============================================================'
write(*,'(A)') ' THERMOFLUIDCALC — SPRAY & DROPLET ATOMIZATION ENGINE'
write(*,'(A)') '============================================================'
write(*,'(A,ES12.4,A)') 'Nozzle Diameter d0 = ', d_nozzle_mm, ' mm'
write(*,'(A,ES12.4,A)') 'Injection Delta-P = ', deltaP_bar, ' bar'
write(*,'(A,ES12.4)') 'Discharge Coefficient Cd = ', Cd
write(*,'(A,ES12.4,A)') 'Liquid Density = ', rho_L, ' kg/m3'
write(*,'(A,ES12.4,A)') 'Liquid Viscosity = ', mu_L, ' Pa.s'
write(*,'(A,ES12.4,A)') 'Surface Tension = ', sigma, ' N/m'
write(*,'(A,ES12.4,A)') 'Gas Density = ', rho_g, ' kg/m3'
write(*,'(A)') '------------------------------------------------------------'
write(*,'(A,ES12.4,A)') 'Injection Velocity = ', V_inj, ' m/s'
write(*,'(A,ES12.4,A)') 'Mass Flow Rate = ', mdot_L, ' kg/s'
write(*,'(A,ES12.4,A)') 'Volumetric Flow Rate = ', Q_L, ' m3/s'
write(*,'(A,ES12.4)') 'Liquid Reynolds Number = ', Re_L
write(*,'(A,ES12.4)') 'Gas Weber Number Weg = ', We_g
write(*,'(A,ES12.4)') 'Ohnesorge Number Oh = ', Oh
write(*,'(A,A)') 'Breakup Regime = ', trim(regime_name)
write(*,'(A)') '------------------------------------------------------------'
write(*,'(A,ES12.4,A)') 'Sauter Mean Diameter D32 = ', D32, ' um'
write(*,'(A,ES12.4,A)') 'Arithmetic Mean D10 = ', D10, ' um'
write(*,'(A,ES12.4,A)') 'Mass Median Diameter Dv50 = ', Dv50, ' um'
write(*,'(A,ES12.4,A)') 'Spray Cone Angle = ', spray_angle, ' deg'
write(*,'(A,ES12.4,A)') 'Liquid Core Breakup Length= ', breakup_length, ' mm'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- ROSIN-RAMMLER DROPLET SIZE DISTRIBUTION ---'
write(*,'(A)') ' d_um Q_cum(%) PDF'
n_points = 25
do i = 0, n_points
d_val = (Dv50 * 2.8d0) * (dble(i) / dble(n_points))
if (d_val <= 1.0d-6) then
cum_vol = 0.0d0
pdf_val = 0.0d0
else
cum_vol = 1.0d0 - exp(-((d_val / X_rr)**q_rr))
pdf_val = (q_rr / X_rr) * ((d_val / X_rr)**(q_rr - 1.0d0)) * exp(-((d_val / X_rr)**q_rr))
end if
write(*,'(2X,F10.2,2X,F10.2,2X,ES12.4)') d_val, cum_vol * 100.0d0, pdf_val * 100.0d0
end do
write(*,*)
write(*,'(A)') '--- SMD D32 VS INJECTION PRESSURE ---'
write(*,'(A)') ' DeltaP_bar D32_um'
do i = 1, 15
p_var = (deltaP_bar * 0.2d0) + (deltaP_bar * 2.8d0) * (dble(i - 1) / 14.0d0)
p_var_Pa = p_var * 1.0d5
m_var = Cd * A_orifice * sqrt(2.0d0 * rho_L * p_var_Pa)
smd_var = (term1 * (m_var**0.22d0) * (p_var_Pa**(-0.43d0))) * 1.0d6
write(*,'(2X,F10.2,2X,F10.2)') p_var, max(1.0d0, smd_var)
end do
end program spray_atomization
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 spray_atomization.f90 -o spray_atomization
2. Execution with input.txt redirection:
spray_atomization < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.50 20.0 0.78 998.0 1.002e-3 0.0728 1.225
Parameter Description:
Nozzle Orifice Diameter [mm] Injection Delta-P [bar] Discharge Coefficient Cd Liquid Density [kg/m3] Liquid Dynamic Viscosity [Pa-s] Surface Tension [N/m] Ambient Gas Density [kg/m3]