🔬
Solver Purpose & Physical Scope
Solves two-phase liquid-gas pressure drops in inclined, horizontal, and vertical pipes using the Friedel and Lockhart-Martinelli two-phase friction multipliers. Computes void fraction and in-situ liquid holdup using Rouhani-Axelsson drift-flux model, separates frictional and hydrostatic gravity head losses, and identifies Taitel-Dukler flow regimes (bubbly, intermittent slug, stratified wavy, annular mist).
📂 Discipline: Fluid
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 390 times
📄 Source File:
two_phase_flow.f90
📁 calcul/FluidMechanics /
two_phase_flow.f90
program two_phase_flow
implicit none
integer :: i, iostat_val, n_points
double precision :: D_mm, L_m, eps_mm, theta_deg, mdot_total, x_quality
double precision :: rho_L, mu_L, sigma, rho_G, mu_G
double precision :: D_m, eps_m, theta_rad, Area, G_flux
double precision :: j_L, j_G, j_tot, Re_L0, Re_G0, f_L0, f_G0
double precision :: rho_H, Fr_H, We_L, E_fr, F_fr, H_fr, phi2_friedel
double precision :: X_tt, phi2_lockhart, alpha_void, HL_holdup, rho_2phi
double precision :: dp_dz_L0, dp_fric_Pa, dp_grav_Pa, dp_total_Pa
double precision :: x_var, x_tt_var, p2_var, void_var
double precision, parameter :: PI = 3.141592653589793d0
double precision, parameter :: g = 9.80665d0
character(len=64) :: regime_name
! Read inputs
read(*,*,iostat=iostat_val) D_mm
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid pipe diameter.'
stop
end if
read(*,*,iostat=iostat_val) L_m
read(*,*,iostat=iostat_val) eps_mm
read(*,*,iostat=iostat_val) theta_deg
read(*,*,iostat=iostat_val) mdot_total
read(*,*,iostat=iostat_val) x_quality
read(*,*,iostat=iostat_val) rho_L
read(*,*,iostat=iostat_val) mu_L
read(*,*,iostat=iostat_val) sigma
read(*,*,iostat=iostat_val) rho_G
read(*,*,iostat=iostat_val) mu_G
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all two-phase flow inputs.'
stop
end if
! Validation
if (D_mm <= 0.0d0 .or. L_m <= 0.0d0 .or. mdot_total <= 0.0d0) then
write(*,*) 'ERROR: Diameter, length, and mass flow rate must be positive.'
stop
end if
if (x_quality < 0.0d0 .or. x_quality > 1.0d0) then
write(*,*) 'ERROR: Vapor quality x must be between 0.0 and 1.0.'
stop
end if
if (rho_L <= 0.0d0 .or. mu_L <= 0.0d0 .or. rho_G <= 0.0d0 .or. mu_G <= 0.0d0) then
write(*,*) 'ERROR: Liquid and gas properties must be positive.'
stop
end if
! Conversions
D_m = D_mm * 1.0d-3
eps_m = eps_mm * 1.0d-3
theta_rad = theta_deg * PI / 180.0d0
Area = (PI / 4.0d0) * (D_m**2)
G_flux = mdot_total / Area ! kg/(m2.s)
! Superficial Velocities
j_L = (G_flux * (1.0d0 - x_quality)) / rho_L
j_G = (G_flux * x_quality) / rho_G
j_tot = j_L + j_G
! Single-phase All-Liquid and All-Gas Reynolds & Friction factors
Re_L0 = (G_flux * D_m) / mu_L
Re_G0 = (G_flux * D_m) / mu_G
! Churchill friction factor for all-liquid
f_L0 = churchill_friction(Re_L0, eps_m / D_m)
f_G0 = churchill_friction(Re_G0, eps_m / D_m)
! Single-Phase Liquid Pressure Gradient (dP/dz)_L0
dp_dz_L0 = f_L0 * (G_flux**2) / (2.0d0 * rho_L * D_m)
! Homogeneous Density & Non-dimensional groups for Friedel
if (x_quality > 0.0d0 .and. x_quality < 1.0d0) then
rho_H = 1.0d0 / ((x_quality / rho_G) + ((1.0d0 - x_quality) / rho_L))
else if (x_quality <= 0.0d0) then
rho_H = rho_L
else
rho_H = rho_G
end if
Fr_H = (G_flux**2) / (g * D_m * (rho_H**2))
if (sigma > 0.0d0) then
We_L = (G_flux**2 * D_m) / (sigma * rho_H)
else
We_L = 100.0d0
end if
! Friedel Two-Phase Multiplier Phi_L0^2
E_fr = ((1.0d0 - x_quality)**2) + (x_quality**2) * ((rho_L * f_G0) / (rho_G * f_L0))
F_fr = (x_quality**0.78d0) * ((1.0d0 - x_quality)**0.224d0)
H_fr = ((rho_L / rho_G)**0.91d0) * ((mu_G / mu_L)**0.19d0) * ((1.0d0 - (mu_G / mu_L))**0.7d0)
if (x_quality > 0.0d0 .and. x_quality < 1.0d0) then
phi2_friedel = E_fr + (3.24d0 * F_fr * H_fr) / ((max(1.0d-4, Fr_H)**0.045d0) * (max(1.0d-4, We_L)**0.035d0))
else
phi2_friedel = 1.0d0
end if
! Lockhart-Martinelli Parameter X_tt
if (x_quality > 0.001d0 .and. x_quality < 0.999d0) then
X_tt = (((1.0d0 - x_quality) / x_quality)**0.9d0) * &
((rho_G / rho_L)**0.5d0) * ((mu_L / mu_G)**0.1d0)
phi2_lockhart = 1.0d0 + (20.0d0 / max(1.0d-3, X_tt)) + (1.0d0 / max(1.0d-6, X_tt**2))
else if (x_quality <= 0.001d0) then
X_tt = 1000.0d0
phi2_lockhart = 1.0d0
else
X_tt = 0.001d0
phi2_lockhart = (rho_L * f_G0) / (rho_G * f_L0)
end if
! Void Fraction (Rouhani-Axelsson Drift Flux Model)
if (x_quality > 0.0d0 .and. x_quality < 1.0d0) then
alpha_void = (x_quality / rho_G) / &
((1.0d0 + 0.12d0 * (1.0d0 - x_quality)) * ((x_quality / rho_G) + ((1.0d0 - x_quality) / rho_L)) + &
(1.18d0 * (1.0d0 - x_quality) * ((g * sigma * (rho_L - rho_G))**0.25d0)) / (G_flux * sqrt(rho_L)))
alpha_void = min(0.999d0, max(0.001d0, alpha_void))
else if (x_quality <= 0.0d0) then
alpha_void = 0.0d0
else
alpha_void = 1.0d0
end if
HL_holdup = 1.0d0 - alpha_void
rho_2phi = (alpha_void * rho_G) + (HL_holdup * rho_L)
! Pressure Drops over Length L
dp_fric_Pa = phi2_friedel * dp_dz_L0 * L_m
dp_grav_Pa = rho_2phi * g * L_m * sin(theta_rad)
dp_total_Pa = dp_fric_Pa + dp_grav_Pa
! Flow Regime Identification
if (alpha_void < 0.20d0 .and. j_L > 0.5d0) then
regime_name = 'Bubbly Flow'
else if (alpha_void < 0.70d0 .and. j_tot < 4.0d0) then
if (abs(theta_deg) < 15.0d0) then
regime_name = 'Stratified / Plug Flow'
else
regime_name = 'Intermittent / Slug Flow'
end if
else if (alpha_void >= 0.70d0 .and. x_quality < 0.85d0) then
regime_name = 'Annular-Mist Flow'
else
regime_name = 'Dispersed Mist Flow'
end if
! Output Formatted Report
write(*,'(A)') '============================================================'
write(*,'(A)') ' THERMOFLUIDCALC — TWO-PHASE PIPE FLOW SOLVER'
write(*,'(A)') '============================================================'
write(*,'(A,ES12.4,A)') 'Pipe Diameter D = ', D_mm, ' mm'
write(*,'(A,ES12.4,A)') 'Pipe Length L = ', L_m, ' m'
write(*,'(A,ES12.4,A)') 'Inclination Angle = ', theta_deg, ' deg'
write(*,'(A,ES12.4,A)') 'Total Mass Flow mdot = ', mdot_total, ' kg/s'
write(*,'(A,ES12.4,A)') 'Mass Flux G = ', G_flux, ' kg/(m2.s)'
write(*,'(A,ES12.4)') 'Vapor/Gas Quality x = ', x_quality
write(*,'(A)') '------------------------------------------------------------'
write(*,'(A,ES12.4,A)') 'Superficial Liquid Vel jL = ', j_L, ' m/s'
write(*,'(A,ES12.4,A)') 'Superficial Gas Vel jG = ', j_G, ' m/s'
write(*,'(A,ES12.4)') 'All-Liquid Reynolds ReL0 = ', Re_L0
write(*,'(A,ES12.4)') 'All-Gas Reynolds ReG0 = ', Re_G0
write(*,'(A,ES12.4)') 'Lockhart-Martinelli Xtt = ', X_tt
write(*,'(A,ES12.4)') 'Friedel Multiplier Phi2L0 = ', phi2_friedel
write(*,'(A)') '------------------------------------------------------------'
write(*,'(A,ES12.4)') 'Void Fraction alpha = ', alpha_void
write(*,'(A,ES12.4)') 'Liquid Holdup HL = ', HL_holdup
write(*,'(A,ES12.4,A)') 'In-situ Mixture Density = ', rho_2phi, ' kg/m3'
write(*,'(A,A)') 'Predicted Flow Regime = ', trim(regime_name)
write(*,'(A)') '------------------------------------------------------------'
write(*,'(A,ES12.4,A)') 'Frictional Pressure Drop = ', dp_fric_Pa * 1.0d-3, ' kPa'
write(*,'(A,ES12.4,A)') 'Hydrostatic Pressure Drop = ', dp_grav_Pa * 1.0d-3, ' kPa'
write(*,'(A,ES12.4,A)') 'Total Two-Phase Delta-P = ', dp_total_Pa * 1.0d-3, ' kPa'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- TWO-PHASE MULTIPLIER AND VOID FRACTION VS QUALITY ---'
write(*,'(A)') ' Quality_x Phi2_Friedel Void_Fraction'
n_points = 20
do i = 1, n_points
x_var = dble(i) / dble(n_points + 1)
E_fr = ((1.0d0 - x_var)**2) + (x_var**2) * ((rho_L * f_G0) / (rho_G * f_L0))
F_fr = (x_var**0.78d0) * ((1.0d0 - x_var)**0.224d0)
p2_var = E_fr + (3.24d0 * F_fr * H_fr) / ((max(1.0d-4, Fr_H)**0.045d0) * (max(1.0d-4, We_L)**0.035d0))
void_var = (x_var / rho_G) / &
((1.0d0 + 0.12d0 * (1.0d0 - x_var)) * ((x_var / rho_G) + ((1.0d0 - x_var) / rho_L)) + &
(1.18d0 * (1.0d0 - x_var) * ((g * sigma * (rho_L - rho_G))**0.25d0)) / (G_flux * sqrt(rho_L)))
write(*,'(2X,F10.3,2X,F12.4,2X,F12.4)') x_var, max(1.0d0, p2_var), min(1.0d0, max(0.0d0, void_var))
end do
contains
function churchill_friction(Re, rel_rough) result(f)
double precision, intent(in) :: Re, rel_rough
double precision :: f, A_term, B_term
if (Re < 1.0d-4) then
f = 0.02d0
return
end if
A_term = (2.457d0 * log(1.0d0 / ((7.0d0 / Re)**0.9d0 + 0.27d0 * rel_rough)))**16
B_term = (37530.0d0 / Re)**16
f = 8.0d0 * (( (8.0d0 / Re)**12 + 1.0d0 / ((A_term + B_term)**1.5d0) )**(1.0d0 / 12.0d0))
end function churchill_friction
end program two_phase_flow
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 two_phase_flow.f90 -o two_phase_flow
2. Execution with input.txt redirection:
two_phase_flow < input.txt
📄 Sample input.txt File Structure
Sample Data:
50.0 25.0 0.045 90.0 2.0 0.20 998.0 1.002e-3 0.0728 1.225 1.81e-5
Parameter Description:
Pipe Diameter D [mm] Pipe Length L [m] Pipe Roughness eps [mm] Inclination Angle theta [deg] (0=horiz, 90=vert) Total Mass Flow Rate [kg/s] Vapor/Gas Quality x [0-1] Liquid Density [kg/m3] Liquid Dynamic Viscosity [Pa-s] Interfacial Surface Tension [N/m] Gas Density [kg/m3] Gas Dynamic Viscosity [Pa-s]