💻 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.
Hardy-Cross Pipe Network Solver
Core Numerical Engine in Fortran 90 • 0 total downloads
pipe_networks.f90
! =========================================================================
! Source File: pipe_networks.f90
! =========================================================================
program pipe_networks
implicit none
integer, parameter :: max_pipes = 8
double precision :: rho, mu, nu, Q_total, dZ
double precision :: L(max_pipes), D(max_pipes), eps(max_pipes)
double precision :: Q(max_pipes), V(max_pipes), Re(max_pipes), f(max_pipes)
double precision :: dP(max_pipes), hf(max_pipes)
double precision :: rel_rough(max_pipes), Ac(max_pipes)
integer :: n_pipes, net_type ! 1=series, 2=parallel
double precision :: pi, g
! Hardy-Cross variables
double precision :: dQ, sum_num, sum_den
double precision :: f_old, Q_old(max_pipes)
integer :: i, j, iter, max_iter, converged
double precision :: tol
! Total results
double precision :: dP_total, hf_total, pump_power
pi = 3.14159265358979d0
g = 9.81d0
max_iter = 500
tol = 1.0d-8
! Read inputs
read(*,*) net_type ! 1=series, 2=parallel
read(*,*) n_pipes ! Number of pipes
read(*,*) rho ! Density [kg/m3]
read(*,*) mu ! Dynamic viscosity [Pa.s]
read(*,*) Q_total ! Total flow rate [m3/s]
read(*,*) dZ ! Elevation difference [m]
nu = mu / rho
do i = 1, n_pipes
read(*,*) L(i), D(i), eps(i)
! D and eps are in meters already (converted by PHP)
Ac(i) = pi * D(i)**2 / 4.0d0
rel_rough(i) = eps(i) / D(i)
end do
! ============================================
! HEADER
! ============================================
write(*,'(A)') '============================================================'
write(*,'(A)') ' PIPE NETWORK ANALYSIS'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUT PARAMETERS ---'
if (net_type == 1) then
write(*,'(A)') ' Network Type = SERIES'
else
write(*,'(A)') ' Network Type = PARALLEL'
end if
write(*,'(A,I2)') ' Number of pipes = ', n_pipes
write(*,'(A,F12.2,A)') ' Fluid density ρ = ', rho, ' kg/m³'
write(*,'(A,ES12.4,A)') ' Dynamic viscosity μ = ', mu, ' Pa·s'
write(*,'(A,ES12.4,A)') ' Kinematic viscosity ν = ', nu, ' m²/s'
write(*,'(A,ES12.4,A)') ' Total flow rate Q = ', Q_total, ' m³/s'
write(*,'(A,F12.4,A)') ' Total flow rate = ', Q_total*1000.0d0*60.0d0, ' L/min'
write(*,'(A,F12.4,A)') ' Elevation difference = ', dZ, ' m'
write(*,*)
write(*,'(A)') '--- PIPE DATA ---'
write(*,'(A)') ' Pipe L [m] D [mm] ε [mm] ε/D A [m²]'
write(*,'(A)') ' ---- ------ ------ ------ ------ ------'
do i = 1, n_pipes
write(*,'(I5,2X,F8.2,2X,F8.2,2X,F8.4,3X,ES10.3,2X,ES10.4)') &
i, L(i), D(i)*1000, eps(i)*1000, rel_rough(i), Ac(i)
end do
write(*,*)
! ============================================
! SERIES NETWORK
! ============================================
if (net_type == 1) then
write(*,'(A)') '--- SERIES ANALYSIS ---'
write(*,'(A)') ' Same flow rate Q in all pipes.'
write(*,*)
dP_total = 0.0d0
hf_total = 0.0d0
do i = 1, n_pipes
Q(i) = Q_total
V(i) = Q(i) / Ac(i)
Re(i) = rho * V(i) * D(i) / mu
! Friction factor
call compute_friction(Re(i), rel_rough(i), f(i))
! Pressure drop (Darcy-Weisbach)
dP(i) = f(i) * (L(i) / D(i)) * rho * V(i)**2 / 2.0d0
hf(i) = dP(i) / (rho * g)
dP_total = dP_total + dP(i)
hf_total = hf_total + hf(i)
end do
! Add elevation head
dP_total = dP_total + rho * g * dZ
hf_total = hf_total + dZ
write(*,'(A)') '--- PER-PIPE RESULTS ---'
write(*,'(A)') ' Pipe Q [L/min] V [m/s] Re f ΔP [Pa] h_f [m]'
write(*,'(A)') ' ---- -------- ------ ---------- -------- --------- -------'
do i = 1, n_pipes
write(*,'(I5,2X,F10.2,2X,F8.4,3X,ES10.4,2X,ES10.4,2X,F10.2,3X,F8.4)') &
i, Q(i)*1000*60, V(i), Re(i), f(i), dP(i), hf(i)
end do
write(*,*)
pump_power = dP_total * Q_total
write(*,'(A)') '--- TOTAL SERIES RESULTS ---'
write(*,'(A,F12.2,A)') ' Total pressure drop = ', dP_total, ' Pa'
write(*,'(A,F12.4,A)') ' Total pressure drop = ', dP_total/1000, ' kPa'
write(*,'(A,F12.4,A)') ' Total head loss = ', hf_total, ' m'
write(*,'(A,F12.4,A)') ' Elevation component = ', rho*g*dZ, ' Pa'
write(*,'(A,F12.4,A)') ' Required pump power = ', pump_power, ' W'
write(*,*)
! ============================================
! PARALLEL NETWORK — HARDY-CROSS METHOD
! ============================================
else
write(*,'(A)') '--- PARALLEL ANALYSIS (HARDY-CROSS) ---'
write(*,'(A)') ' Same pressure drop across all branches.'
write(*,'(A)') ' Iterating to find flow distribution...'
write(*,*)
! Initial guess: distribute proportional to D^2.5
sum_num = 0.0d0
do i = 1, n_pipes
sum_num = sum_num + D(i)**2.5d0 / L(i)**0.5d0
end do
do i = 1, n_pipes
Q(i) = Q_total * (D(i)**2.5d0 / L(i)**0.5d0) / sum_num
end do
write(*,'(A)') '--- CONVERGENCE HISTORY ---'
write(*,'(A)') ' Iter ΔQ [m³/s] Max |Q change| Σ Q [m³/s]'
write(*,'(A)') ' ---- ---------- --------------- ----------'
converged = 0
do iter = 1, max_iter
! Save old values
do i = 1, n_pipes
Q_old(i) = Q(i)
end do
! Compute friction factors and head losses
do i = 1, n_pipes
V(i) = Q(i) / Ac(i)
Re(i) = rho * abs(V(i)) * D(i) / mu
if (Re(i) < 1.0d0) Re(i) = 1.0d0
call compute_friction(Re(i), rel_rough(i), f(i))
end do
! Hardy-Cross correction for parallel network
! For equal head loss: h_f1 = h_f2 = ... = h_fn
! Use head loss formulation: h_f = (8 f L Q|Q|) / (π² g D⁵)
! Correction: ΔQ = -Σ(R_i Q_i |Q_i|) / Σ(2 R_i |Q_i|)
! where R_i = 8 f_i L_i / (π² g D_i⁵)
! Actually for parallel, we equalize head losses
! Approach: adjust Q's so all have same head loss
! Using the Newton-like correction:
! Compute head losses
do i = 1, n_pipes
hf(i) = f(i) * L(i) * V(i) * abs(V(i)) / (2.0d0 * g * D(i))
end do
! Target: average head loss
sum_num = 0.0d0
do i = 1, n_pipes
sum_num = sum_num + hf(i)
end do
sum_num = sum_num / dble(n_pipes) ! average h_f
! Adjust each pipe's flow to equalize head loss
! h_f = f L V² / (2gD) = f L Q² / (2g D A²)
! dh_f/dQ = f L 2Q / (2g D A²) = 2 h_f / Q (approximately)
! ΔQ_i = -(h_f_i - h_f_avg) / (2 h_f_i / Q_i)
do i = 1, n_pipes
if (abs(hf(i)) > 1.0d-15 .and. abs(Q(i)) > 1.0d-15) then
dQ = -(hf(i) - sum_num) / (2.0d0 * hf(i) / Q(i))
Q(i) = Q(i) + dQ * 0.5d0 ! relaxation factor
end if
end do
! Enforce mass conservation: Σ Q_i = Q_total
sum_den = 0.0d0
do i = 1, n_pipes
sum_den = sum_den + Q(i)
end do
do i = 1, n_pipes
Q(i) = Q(i) * Q_total / sum_den
end do
! Check convergence
dQ = 0.0d0
do i = 1, n_pipes
if (abs(Q(i) - Q_old(i)) > dQ) dQ = abs(Q(i) - Q_old(i))
end do
if (mod(iter, 10) == 1 .or. iter <= 5 .or. dQ < tol) then
sum_den = 0.0d0
do i = 1, n_pipes
sum_den = sum_den + Q(i)
end do
write(*,'(I5,3X,ES14.6,3X,ES14.6,3X,ES12.4)') iter, dQ, dQ, sum_den
end if
if (dQ < tol) then
converged = 1
write(*,*)
write(*,'(A,I4,A)') ' ✓ Converged after ', iter, ' iterations'
exit
end if
end do
if (converged == 0) then
write(*,*)
write(*,'(A)') ' ⚠ WARNING: Did not converge within max iterations'
end if
write(*,*)
! Final results
do i = 1, n_pipes
V(i) = Q(i) / Ac(i)
Re(i) = rho * abs(V(i)) * D(i) / mu
call compute_friction(Re(i), rel_rough(i), f(i))
dP(i) = f(i) * (L(i) / D(i)) * rho * V(i)**2 / 2.0d0
hf(i) = dP(i) / (rho * g)
end do
write(*,'(A)') '--- PER-PIPE RESULTS ---'
write(*,'(A)') ' Pipe Q [L/min] V [m/s] Re f ΔP [Pa] h_f [m] Q/Q_t [%]'
write(*,'(A)') ' ---- -------- ------ ---------- -------- --------- ------- ---------'
do i = 1, n_pipes
write(*,'(I5,2X,F10.2,2X,F8.4,3X,ES10.4,2X,ES10.4,2X,F10.2,3X,F8.4,3X,F8.2)') &
i, Q(i)*1000*60, V(i), Re(i), f(i), dP(i), hf(i), (Q(i)/Q_total)*100
end do
write(*,*)
! Use first pipe's head loss as reference (all should be equal)
dP_total = dP(1)
hf_total = hf(1)
pump_power = dP_total * Q_total
write(*,'(A)') '--- TOTAL PARALLEL RESULTS ---'
write(*,'(A,F12.2,A)') ' Common pressure drop = ', dP_total, ' Pa'
write(*,'(A,F12.4,A)') ' Common pressure drop = ', dP_total/1000, ' kPa'
write(*,'(A,F12.4,A)') ' Common head loss = ', hf_total, ' m'
write(*,'(A,F12.4,A)') ' Required pump power = ', pump_power, ' W'
write(*,'(A,ES12.4,A)') ' Total flow rate check = ', Q_total, ' m³/s'
write(*,*)
end if
! ============================================
! CHART DATA
! ============================================
write(*,'(A)') '--- CHART_DATA ---'
do i = 1, n_pipes
write(*,'(I2,A,F12.6,A,F12.4,A,F12.2,A,F12.6)') &
i, ',', Q(i)*1000*60, ',', V(i), ',', dP(i), ',', hf(i)
end do
write(*,'(A)') '--- END_CHART_DATA ---'
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ---'
write(*,'(A)') ' Darcy-Weisbach: ΔP = f(L/D)(ρV²/2)'
write(*,'(A)') ' Colebrook-White: 1/√f = -2log₁₀(ε/D/3.7 + 2.51/(Re√f))'
write(*,'(A)') ' Series: Q₁ = Q₂ = ... = Q_n; ΔP_total = Σ ΔP_i'
write(*,'(A)') ' Parallel: ΔP₁ = ΔP₂ = ... = ΔP_n; Q_total = Σ Q_i'
write(*,'(A)') ' Hardy-Cross: ΔQ = -Σ(RQ|Q|) / Σ(2R|Q|)'
contains
subroutine compute_friction(Re_val, eD, f_out)
double precision, intent(in) :: Re_val, eD
double precision, intent(out) :: f_out
double precision :: f_old_loc
integer :: k
if (Re_val < 2300.0d0) then
f_out = 64.0d0 / Re_val
else
! Swamee-Jain initial guess
f_out = 0.25d0 / (log10(eD/3.7d0 + 5.74d0/Re_val**0.9d0))**2
! Colebrook-White iteration
do k = 1, 100
f_old_loc = f_out
f_out = 1.0d0 / (-2.0d0 * log10(eD/3.7d0 + 2.51d0/(Re_val*sqrt(f_old_loc))))**2
if (abs(f_out - f_old_loc) / f_out < 1.0d-10) exit
end do
end if
end subroutine
end program pipe_networks
Solver Description
Solves steady pipe networks using the Hardy-Cross loop-balancing method. Iteratively solves for flow rates in loops by setting $\Delta Q = -\sum h_f / (2 \sum (h_f / Q))$ until mass conservation at nodes and energy conservation in loops are satisfied.
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 pipe_networks.f90 -o pipe_networks_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
pipe_networks_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Fluid Density ($\rho$) [kg/m³]
1000.0
! Fluid Viscosity ($\mu$) [Pa-s]
0.001
! Number of loops (N_loops)
2
! Boundary Flow Rate [L/min]
10.0
! Elevation Difference $\Delta Z$ [m]
0.0
! Pipe Index 1
1
! Pipe 1 Length [m]
100.0
! Pipe 1 Diameter [mm]
100.0
! Pipe 1 Roughness [mm]
0.045
! Pipe Index 2
2
! Pipe 2 Length [m]
80.0
! Pipe 2 Diameter [mm]
80.0
! Pipe 2 Roughness [mm]
0.15
1000.0
! Fluid Viscosity ($\mu$) [Pa-s]
0.001
! Number of loops (N_loops)
2
! Boundary Flow Rate [L/min]
10.0
! Elevation Difference $\Delta Z$ [m]
0.0
! Pipe Index 1
1
! Pipe 1 Length [m]
100.0
! Pipe 1 Diameter [mm]
100.0
! Pipe 1 Roughness [mm]
0.045
! Pipe Index 2
2
! Pipe 2 Length [m]
80.0
! Pipe 2 Diameter [mm]
80.0
! Pipe 2 Roughness [mm]
0.15