💻 Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Pipe Network Analysis
Core Numerical Engine in Fortran 90 • 85 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
Pipe networks are solved using the conservation of mass at nodes and conservation of energy in loops (sum of head losses is zero). The Hardy-Cross method iteratively computes flow rate corrections $\Delta Q$ for each loop:
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 -O3 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):
! Network type (1=Loop, 2=Branch)
1
! Number of pipes
3
! Fluid density [kg/m3]
998.0
! Fluid dynamic viscosity [Pa-s]
0.001002
! Total inlet flow rate [m3/s]
0.01
! Elevation change [m]
0.0
! Pipe definitions (Loop_Node1 Loop_Node2 Diam_m Len_m Rough_m K_minor)
1 2 0.05 100.0 4.5e-5 1.5
! Parameter 8
2 3 0.04 150.0 4.5e-5 2.0
! Parameter 9
3 1 0.05 120.0 4.5e-5 0.5
1
! Number of pipes
3
! Fluid density [kg/m3]
998.0
! Fluid dynamic viscosity [Pa-s]
0.001002
! Total inlet flow rate [m3/s]
0.01
! Elevation change [m]
0.0
! Pipe definitions (Loop_Node1 Loop_Node2 Diam_m Len_m Rough_m K_minor)
1 2 0.05 100.0 4.5e-5 1.5
! Parameter 8
2 3 0.04 150.0 4.5e-5 2.0
! Parameter 9
3 1 0.05 120.0 4.5e-5 0.5