🔬
Solver Purpose & Physical Scope
Flow in pipes experiences pressure drops due to viscous friction and minor losses (valves, bends). Friction factor is determined using the Colebrook-White equation for turbulent flow ($Re > 4000$):
📂 Discipline: Fluid
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 302 times
📄 Source File:
pipe_flow_calculator.f90
📁 calcul/FluidMechanics /
pipe_flow_calculator.f90
program pipe_flow_calculator
implicit none
double precision :: D, L, roughness, rho, mu, nu
double precision :: Q_flow, U_m, Re, f, f_old, dP, dP_kPa
double precision :: A_c, pi, V_dot
double precision :: head_loss, pump_power
double precision :: f_lam, f_turb_smooth, f_turb_rough
double precision :: rel_rough
integer :: i, n_iter, converged
character(len=30) :: flow_regime
! Minor loss variables
integer :: n_fittings
double precision :: K_total, dP_minor, head_minor
pi = 3.14159265358979d0
! Read inputs
read(*,*) D ! Pipe inner diameter [m]
read(*,*) L ! Pipe length [m]
read(*,*) roughness ! Surface roughness eps [m]
read(*,*) rho ! Fluid density [kg/m3]
read(*,*) mu ! Dynamic viscosity [Pa.s]
read(*,*) V_dot ! Volume flow rate [m3/s]
read(*,*) K_total ! Sum of minor loss coefficients
! Computed properties
nu = mu / rho
A_c = pi * D**2 / 4.0d0
U_m = V_dot / A_c
Re = rho * U_m * D / mu
rel_rough = roughness / D
! ============================================
! FLOW REGIME DETERMINATION
! ============================================
if (Re < 2300.0d0) then
flow_regime = 'Laminar'
else if (Re < 4000.0d0) then
flow_regime = 'Transitional'
else
flow_regime = 'Turbulent'
end if
! ============================================
! FRICTION FACTOR CALCULATION
! ============================================
! Laminar friction factor
f_lam = 64.0d0 / Re
! Turbulent: Colebrook-White equation (iterative)
! 1/√f = -2.0 log1₀(eps/D/3.7 + 2.51/(Re√f))
if (Re >= 2300.0d0) then
! Initial guess: Swamee-Jain approximation
f = 0.25d0 / (log10(rel_rough/3.7d0 + 5.74d0/Re**0.9d0))**2
converged = 0
n_iter = 0
do i = 1, 100
f_old = f
f = 1.0d0 / (-2.0d0 * log10(rel_rough/3.7d0 + 2.51d0/(Re*sqrt(f_old))))**2
n_iter = i
if (abs(f - f_old) / f < 1.0d-8) then
converged = 1
exit
end if
end do
f_turb_rough = f
else
f = f_lam
f_turb_rough = 0.0d0
n_iter = 0
converged = 1
end if
! Smooth pipe (Petukhov)
if (Re >= 2300.0d0) then
f_turb_smooth = (0.790d0 * log(Re) - 1.64d0)**(-2.0d0)
else
f_turb_smooth = 0.0d0
end if
! ============================================
! PRESSURE DROP (Darcy-Weisbach)
! ============================================
! deltaP = f x (L/D) x (rhoU2/2)
dP = f * (L / D) * rho * U_m**2 / 2.0d0
dP_kPa = dP / 1000.0d0
head_loss = dP / (rho * 9.81d0)
! Minor losses
dP_minor = K_total * rho * U_m**2 / 2.0d0
head_minor = dP_minor / (rho * 9.81d0)
! Total
pump_power = (dP + dP_minor) * V_dot
! ============================================
! OUTPUT
! ============================================
write(*,'(A)') '============================================================'
write(*,'(A)') ' PIPE FLOW CALCULATOR — DARCY-WEISBACH METHOD'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- PIPE GEOMETRY ---'
write(*,'(A,F12.2,A)') ' Inner Diameter (D) = ', D*1000, ' mm'
write(*,'(A,F12.4,A)') ' Pipe Length (L) = ', L, ' m'
write(*,'(A,ES12.4,A)') ' Surface Roughness (eps) = ', roughness, ' m'
write(*,'(A,ES12.4)') ' Relative Roughness eps/D = ', rel_rough
write(*,'(A,ES12.4,A)') ' Cross-Section Area = ', A_c, ' m2'
write(*,*)
write(*,'(A)') '--- FLUID PROPERTIES ---'
write(*,'(A,F12.2,A)') ' Density (rho) = ', rho, ' kg/m3'
write(*,'(A,ES12.4,A)') ' Dynamic Viscosity (mu) = ', mu, ' Pa.s'
write(*,'(A,ES12.4,A)') ' Kinematic Viscosity (nu) = ', nu, ' m2/s'
write(*,*)
write(*,'(A)') '--- FLOW CONDITIONS ---'
write(*,'(A,ES12.4,A)') ' Volume Flow Rate (Q) = ', V_dot, ' m3/s'
write(*,'(A,F12.4,A)') ' Volume Flow Rate = ', V_dot*1000*60, ' L/min'
write(*,'(A,F12.4,A)') ' Mean Velocity (U) = ', U_m, ' m/s'
write(*,'(A,ES12.4)') ' Reynolds Number (Re) = ', Re
write(*,'(A,A)') ' Flow Regime = ', trim(flow_regime)
write(*,*)
write(*,'(A)') '--- FRICTION FACTOR ---'
write(*,'(A,ES12.6)') ' f (Darcy) = ', f
write(*,'(A,ES12.6)') ' f_laminar (64/Re) = ', f_lam
if (Re >= 2300.0d0) then
write(*,'(A,ES12.6)') ' f_smooth (Petukhov) = ', f_turb_smooth
write(*,'(A,ES12.6)') ' f_rough (Colebrook) = ', f_turb_rough
write(*,'(A,I4,A)') ' Colebrook iterations = ', n_iter, ' (converged)'
end if
write(*,*)
write(*,'(A)') '--- MAJOR LOSSES (friction) ---'
write(*,'(A,F12.2,A)') ' Pressure Drop (deltaP) = ', dP, ' Pa'
write(*,'(A,F12.4,A)') ' Pressure Drop = ', dP_kPa, ' kPa'
write(*,'(A,F12.4,A)') ' Head Loss (h_f) = ', head_loss, ' m'
write(*,*)
write(*,'(A)') '--- MINOR LOSSES (fittings) ---'
write(*,'(A,F12.4)') ' Total K coefficient = ', K_total
write(*,'(A,F12.2,A)') ' Minor Pressure Drop = ', dP_minor, ' Pa'
write(*,'(A,F12.4,A)') ' Minor Head Loss = ', head_minor, ' m'
write(*,*)
write(*,'(A)') '--- TOTAL RESULTS ---'
write(*,'(A,F12.2,A)') ' Total Pressure Drop = ', dP + dP_minor, ' Pa'
write(*,'(A,F12.4,A)') ' Total Pressure Drop = ', (dP + dP_minor)/1000, ' kPa'
write(*,'(A,F12.4,A)') ' Total Head Loss = ', head_loss + head_minor, ' m'
write(*,'(A,F12.4,A)') ' Required Pump Power = ', pump_power, ' W'
write(*,*)
! Moody chart data points
write(*,'(A)') '--- FRICTION FACTOR vs Re (for Moody chart) ---'
write(*,'(A)') ' Re f_smooth f_rough f_laminar'
write(*,'(A)') ' ---'
do i = 1, 30
! Logarithmic spacing from Re=500 to Re=1e7
Re = 10.0d0**(2.7d0 + dble(i-1) * 4.3d0 / 29.0d0)
f_lam = 64.0d0 / Re
! Colebrook for rough pipe
f = 0.25d0 / (log10(rel_rough/3.7d0 + 5.74d0/Re**0.9d0))**2
do n_iter = 1, 50
f_old = f
f = 1.0d0 / (-2.0d0 * log10(rel_rough/3.7d0 + 2.51d0/(Re*sqrt(f_old))))**2
if (abs(f - f_old) / f < 1.0d-8) exit
end do
f_turb_rough = f
! Smooth
f_turb_smooth = (0.790d0 * log(Re) - 1.64d0)**(-2.0d0)
if (Re < 2300.0d0) then
write(*,'(ES12.4,2X,F12.8,2X,F12.8,2X,F12.8)') Re, f_turb_smooth, f_turb_rough, f_lam
else
write(*,'(ES12.4,2X,F12.8,2X,F12.8,2X,A)') Re, f_turb_smooth, f_turb_rough, ' ---'
end if
end do
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ---'
write(*,'(A)') ' Darcy-Weisbach: deltaP = f(L/D)(rhoU2/2)'
write(*,'(A)') ' Colebrook: 1/√f = -2log1₀(eps/D/3.7 + 2.51/(Re√f))'
write(*,'(A)') ' Minor losses: deltaP_minor = ΣK(rhoU2/2)'
write(*,'(A)') ' Head loss: h_f = deltaP/(rhog)'
end program pipe_flow_calculator
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 pipe_flow_calculator.f90 -o pipe_flow_calc
2. Execution with input.txt redirection:
pipe_flow_calc < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.05 100.0 4.5e-5 998.0 0.001002 0.005 2.5
Parameter Description:
Pipe inner diameter [m] Pipe length [m] Roughness epsilon [m] Fluid density [kg/m3] Fluid dynamic viscosity [Pa-s] Volumetric flow rate [m3/s] Sum of minor loss coefficients K