Python Project: Pipe Network Solver

In this capstone project we'll combine everything from the course — fluid properties, the Reynolds number, friction factors, and the Darcy-Weisbach equation — into a complete, reusable Python tool. We'll compute head losses for a single pipe, size a pump, and then tackle the real engineering challenge: solving a looped pipe network with the classic Hardy Cross method.

This is exactly the kind of calculation that sits inside commercial software like EPANET, used to design municipal water distribution systems.

The Physics We Are Coding

$$Re = \frac{\rho V D}{\mu} = \frac{V D}{\nu}$$

$$h_f = f \frac{L}{D}\frac{V^2}{2g}$$

The friction factor $f$ comes from the flow regime:

• Laminar ($Re < 2300$): $f = \dfrac{64}{Re}$
• Turbulent: the Colebrook equation (implicit):

$$\frac{1}{\sqrt{f}} = -2\log_{10}\!\left(\frac{\varepsilon/D}{3.7} + \frac{2.51}{Re\sqrt{f}}\right)$$

or the explicit Swamee-Jain approximation:

$$f = \frac{0.25}{\left[\log_{10}\!\left(\dfrac{\varepsilon/D}{3.7} + \dfrac{5.74}{Re^{0.9}}\right)\right]^2}$$

$$h_m = K \frac{V^2}{2g}$$

Project Structure

pipe_network/
├── fluids.py       # Fluid property definitions
├── pipe.py         # Single Pipe class (Re, f, head loss)
├── pump.py         # Pump head sizing
└── network.py      # Hardy Cross looped-network solver

Complete Implementation

import numpy as np
import matplotlib.pyplot as plt
from dataclasses import dataclass, field
from scipy.optimize import brentq

g = 9.81   # gravitational acceleration (m/s^2)

# ============================================================
# FLUIDS MODULE
# ============================================================

@dataclass
class Fluid:
    """Fluid properties at a given temperature."""
    name: str
    rho: float     # density (kg/m^3)
    mu: float      # dynamic viscosity (Pa.s)

    @property
    def nu(self):
        """Kinematic viscosity (m^2/s)."""
        return self.mu / self.rho

    @classmethod
    def water(cls):
        return cls("Water (20C)", rho=998.0, mu=1.002e-3)

    @classmethod
    def oil(cls):
        return cls("Light oil", rho=870.0, mu=0.08)

    @classmethod
    def air(cls):
        return cls("Air (20C)", rho=1.204, mu=1.81e-5)


# ============================================================
# FRICTION FACTOR MODELS
# ============================================================

def colebrook(Re, rel_rough):
    """Solve the implicit Colebrook equation for friction factor f."""
    if Re < 2300:
        return 64.0 / Re
    def residual(f):
        return 1/np.sqrt(f) + 2*np.log10(rel_rough/3.7 + 2.51/(Re*np.sqrt(f)))
    return brentq(residual, 1e-4, 1.0)

def swamee_jain(Re, rel_rough):
    """Explicit Swamee-Jain approximation for turbulent friction factor."""
    if Re < 2300:
        return 64.0 / Re
    return 0.25 / (np.log10(rel_rough/3.7 + 5.74/Re**0.9))**2


# ============================================================
# PIPE MODULE
# ============================================================

@dataclass
class Pipe:
    """A single pipe segment with major and minor losses."""
    length: float            # m
    diameter: float          # m
    roughness: float = 1.5e-6  # m (drawn tubing default)
    K_total: float = 0.0     # sum of minor-loss coefficients
    fluid: Fluid = field(default_factory=Fluid.water)

    @property
    def area(self):
        return np.pi * self.diameter**2 / 4

    def velocity(self, Q):
        """Mean velocity (m/s) for volumetric flow Q (m^3/s)."""
        return Q / self.area

    def reynolds(self, Q):
        V = self.velocity(Q)
        return V * self.diameter / self.fluid.nu

    def friction_factor(self, Q, method="colebrook"):
        Re = self.reynolds(Q)
        rel = self.roughness / self.diameter
        if method == "swamee_jain":
            return swamee_jain(Re, rel)
        return colebrook(Re, rel)

    def head_loss(self, Q, method="colebrook"):
        """Total head loss h_f + h_m (m) for flow Q (m^3/s)."""
        V = self.velocity(Q)
        f = self.friction_factor(abs(Q), method)
        hf = f * (self.length / self.diameter) * V**2 / (2 * g)
        hm = self.K_total * V**2 / (2 * g)
        return hf + hm

    def report(self, Q):
        Re = self.reynolds(Q)
        f = self.friction_factor(Q)
        regime = "laminar" if Re < 2300 else "turbulent"
        print(f"  Q={Q*1000:6.2f} L/s | V={self.velocity(Q):5.2f} m/s | "
              f"Re={Re:9.0f} ({regime}) | f={f:.4f} | "
              f"h_loss={self.head_loss(Q):6.2f} m")


# ============================================================
# PUMP MODULE
# ============================================================

def required_pump_head(pipe, Q, z_static):
    """
    Required pump head (m) = static lift + friction/minor losses.
    z_static : elevation gain from source to delivery (m).
    """
    return z_static + pipe.head_loss(Q)

def pump_power(pipe, Q, z_static, efficiency=0.75):
    """Shaft power (W) for a pump delivering Q against head H."""
    H = required_pump_head(pipe, Q, z_static)
    return pipe.fluid.rho * g * Q * H / efficiency


# ============================================================
# HARDY CROSS NETWORK SOLVER
# ============================================================

class Network:
    """
    Looped pipe network solved by the Hardy Cross method.

    Head loss is modelled as h = r * Q^n with n = 2 (turbulent).
    Each loop is a list of (pipe_index, sign) where sign = +1 if the
    assumed flow direction agrees with the clockwise loop traversal.
    """

    def __init__(self, n=2.0):
        self.pipes = []        # list of dicts: {'r': resistance, 'Q': flow}
        self.loops = []        # list of [(pipe_idx, sign), ...]
        self.n = n

    def add_pipe(self, r, Q0):
        """Add a pipe with resistance r (so h = r*Q^n) and initial guess Q0."""
        self.pipes.append({'r': r, 'Q': Q0})
        return len(self.pipes) - 1

    def add_loop(self, members):
        self.loops.append(members)

    def solve(self, tol=1e-6, max_iter=100):
        for it in range(max_iter):
            max_dQ = 0.0
            for loop in self.loops:
                num = 0.0   # sum of signed head losses
                den = 0.0   # sum of n*|h|/Q
                for idx, sign in loop:
                    Q = self.pipes[idx]['Q']
                    r = self.pipes[idx]['r']
                    h = r * abs(Q)**(self.n - 1) * Q   # signed head loss
                    num += sign * h
                    den += self.n * r * abs(Q)**(self.n - 1)
                dQ = -num / den if den != 0 else 0.0
                # Apply correction to every pipe in the loop
                for idx, sign in loop:
                    self.pipes[idx]['Q'] += sign * dQ
                max_dQ = max(max_dQ, abs(dQ))
            if max_dQ < tol:
                return it + 1
        return max_iter

    def flows(self):
        return [p['Q'] for p in self.pipes]


# ============================================================
# EXAMPLE USAGE
# ============================================================

if __name__ == "__main__":
    print("=" * 64)
    print("SINGLE PIPE ANALYSIS")
    print("=" * 64)
    pipe = Pipe(length=500, diameter=0.30, roughness=0.26e-3,
                K_total=0.5 + 4*0.9 + 1.0, fluid=Fluid.water())
    for Q in [0.05, 0.10, 0.20, 0.30]:   # m^3/s
        pipe.report(Q)

    print("\nPUMP SIZING (deliver 0.20 m^3/s, lift 25 m):")
    Q = 0.20
    H = required_pump_head(pipe, Q, z_static=25)
    P = pump_power(pipe, Q, z_static=25, efficiency=0.75)
    print(f"  Required head  = {H:.2f} m")
    print(f"  Shaft power    = {P/1000:.2f} kW")

================================================================
SINGLE PIPE ANALYSIS
================================================================
  Q=  50.00 L/s | V= 0.71 m/s | Re=   211359 (turbulent) | f=0.0204 | h_loss=  1.00 m
  Q= 100.00 L/s | V= 1.41 m/s | Re=   422719 (turbulent) | f=0.0197 | h_loss=  3.88 m
  Q= 200.00 L/s | V= 2.83 m/s | Re=   845438 (turbulent) | f=0.0194 | h_loss= 15.25 m
  Q= 300.00 L/s | V= 4.24 m/s | Re=  1268157 (turbulent) | f=0.0192 | h_loss= 34.12 m

PUMP SIZING (deliver 0.20 m^3/s, lift 25 m):
  Required head  = 40.25 m
  Shaft power    = 105.09 kW

Comparing Friction Factor Models

The implicit Colebrook equation is the standard, but Swamee-Jain is an explicit one-line approximation. Let's check how close they agree across the turbulent range.

import numpy as np

pipe = Pipe(length=100, diameter=0.20, roughness=0.046e-3, fluid=Fluid.water())

print(f"{'Q (L/s)':>8} {'Re':>10} {'Colebrook':>12} {'Swamee-Jain':>12} {'% diff':>8}")
for Q in [0.02, 0.05, 0.10, 0.20, 0.40]:
    Re = pipe.reynolds(Q)
    fc = pipe.friction_factor(Q, "colebrook")
    fs = pipe.friction_factor(Q, "swamee_jain")
    diff = 100 * (fs - fc) / fc
    print(f"{Q*1000:8.1f} {Re:10.0f} {fc:12.5f} {fs:12.5f} {diff:7.2f}%")

 Q (L/s)         Re    Colebrook  Swamee-Jain   % diff
    20.0     126816      0.01844      0.01846    0.09%
    50.0     317039      0.01638      0.01645    0.43%
   100.0     634078      0.01542      0.01551    0.58%
   200.0    1268157      0.01483      0.01492    0.59%
   400.0    2536314      0.01450      0.01457    0.49%

Swamee-Jain stays within about 0.6% of Colebrook — excellent for hand calculations and fast network solvers.

Solving a Looped Network with Hardy Cross

Consider a single-loop network: an inflow of 100 L/s enters at node A, splits into two parallel paths, and rejoins at node C before leaving (100 L/s out). The clockwise path A→B→C and the counter-clockwise path A→D→C form one loop. We assign each pipe a resistance $r$ (so $h = rQ^2$) and an initial guess that satisfies continuity at every node, then iterate corrections until the head loss around the loop sums to zero.

import numpy as np

# Inflow 100 L/s at A splits into two parallel paths to C.
# Continuity: Q_AB + Q_AD = 0.10 m^3/s  (initial guess 0.06 + 0.04)
net = Network(n=2.0)

# Path 1 (clockwise, sign +1):  A -> B -> C
p_AB = net.add_pipe(r=80.0, Q0=0.06)   # 0
p_BC = net.add_pipe(r=30.0, Q0=0.06)   # 1
# Path 2 (counter-clockwise, sign -1): A -> D -> C
p_AD = net.add_pipe(r=50.0, Q0=0.04)   # 2
p_DC = net.add_pipe(r=20.0, Q0=0.04)   # 3

# One loop: go +path1, -path2 (so head losses balance around the loop)
net.add_loop([(p_AB, +1), (p_BC, +1), (p_AD, -1), (p_DC, -1)])

iters = net.solve()
print(f"Converged in {iters} iterations\n")

names = ["A-B", "B-C", "A-D", "D-C"]
flows = net.flows()
for name, Q in zip(names, flows):
    print(f"  Pipe {name}:  Q = {Q*1000:6.2f} L/s")

# Continuity check at the split: both paths must sum to the 100 L/s inflow
print(f"\nPath A-B-C = {flows[0]*1000:.2f} L/s,  Path A-D-C = {flows[2]*1000:.2f} L/s")
print(f"Total to C = {(flows[0] + flows[2])*1000:.2f} L/s  (inflow = 100.00 L/s)")

# Loop head-loss closure should be ~0
closure = sum(s * net.pipes[i]['r'] * abs(flows[i])*flows[i]
              for i, s in [(0, 1), (1, 1), (2, -1), (3, -1)])
print(f"Loop head-loss closure = {closure:.2e} m  (should be ~0)")

Converged in 4 iterations

  Pipe A-B:  Q =  44.37 L/s
  Pipe B-C:  Q =  44.37 L/s
  Pipe A-D:  Q =  55.63 L/s
  Pipe D-C:  Q =  55.63 L/s

Path A-B-C = 44.37 L/s,  Path A-D-C = 55.63 L/s
Total to C = 100.00 L/s  (inflow = 100.00 L/s)
Loop head-loss closure = 6.94e-18 m  (should be ~0)

The flow self-distributes so the head loss is identical along both paths (closure ~$10^{-18}$) while continuity is preserved — the lower-resistance path D carries more flow. This is Kirchhoff's voltage and current laws applied to pipes.

Head-Loss Visualization

Finally, let's visualize how head loss grows with flow rate for several pipe diameters — the chart an engineer uses to pick a pipe size for a target flow.

import numpy as np
import matplotlib.pyplot as plt

water = Fluid.water()
Q = np.linspace(0.005, 0.30, 120)       # m^3/s
diameters = [0.15, 0.20, 0.25, 0.30]    # m

plt.figure(figsize=(9, 6))
for D in diameters:
    pipe = Pipe(length=500, diameter=D, roughness=0.26e-3,
                K_total=2.0, fluid=water)
    h = np.array([pipe.head_loss(q) for q in Q])
    plt.plot(Q*1000, h, lw=2, label=f'D = {D*1000:.0f} mm')

plt.axhline(40, color='gray', ls='--', alpha=0.7, label='Pump limit (40 m)')
plt.xlabel('Flow rate Q (L/s)')
plt.ylabel('Total head loss (m)')
plt.title('Head Loss vs Flow Rate for a 500 m Pipe')
plt.legend(); plt.grid(alpha=0.3)
plt.ylim(0, 80)
plt.show()

# Largest flow each pipe can carry under a 40 m head budget
print("Max flow within 40 m head budget:")
for D in diameters:
    pipe = Pipe(length=500, diameter=D, roughness=0.26e-3,
                K_total=2.0, fluid=water)
    h = np.array([pipe.head_loss(q) for q in Q])
    q_max = Q[h <= 40][-1] if np.any(h <= 40) else 0
    print(f"  D = {D*1000:.0f} mm  ->  Q_max = {q_max*1000:.1f} L/s")

Max flow within 40 m head budget:
  D = 150 mm  ->  Q_max = 54.6 L/s
  D = 200 mm  ->  Q_max = 116.6 L/s
  D = 250 mm  ->  Q_max = 210.8 L/s
  D = 300 mm  ->  Q_max = 300.0 L/s

Doubling the diameter slashes head loss roughly thirty-fold (head loss $\propto 1/D^5$ for fixed flow) — which is why upsizing a pipe is often cheaper than upsizing a pump.

Extending the Tool

Ideas to take this project further:

• Node-based solver: Solve for nodal heads directly (gradient/Newton method, as EPANET does).
• Pump curves: Replace the fixed pump head with a real $H$-$Q$ characteristic curve.
• Temperature-dependent properties: Interpolate $\rho$ and $\mu$ from data tables.
• Series/parallel reductions: Auto-combine pipes before solving.
• Cost optimization: Minimize pipe + pumping cost over a design lifetime.

Common Pitfalls

• Guessing flows that violate continuity. Hardy Cross requires the initial flows to balance at every node, or it will not converge to a physical answer.
• Wrong loop sign convention. A pipe shared between two loops takes opposite signs in each; mixing this up gives nonsense.
• Using $f = 64/Re$ in turbulent flow. That formula is laminar-only; always branch on $Re = 2300$.
• Ignoring minor losses. For short pipes with many fittings, $h_m$ can exceed $h_f$.
• Forgetting static lift in pump sizing. Required head = elevation gain plus friction losses, not just friction.

Key Takeaways

• A single Pipe class cleanly bundles velocity, Reynolds number, friction factor, and Darcy-Weisbach head loss.
• Colebrook (implicit) and Swamee-Jain (explicit) agree within ~0.3%; use whichever fits your speed/accuracy needs.
• Pump head = static lift + total head loss; power follows from $P = \rho g Q H / \eta$.
• The Hardy Cross method solves looped networks by iteratively correcting loop flows until head losses close to zero.
• Head loss scales steeply with diameter ($\propto 1/D^5$), so pipe sizing is one of the most powerful design levers.

Congratulations — you've built a real engineering tool and completed the Fluid Mechanics with Python course! Take the final quiz to test everything you've learned.