Hardy Cross method of solution of the pipe network

hardy_cross_example.ipynb

{
 "metadata": {
  "name": "hardy_cross_example"
 }, 
 "nbformat": 2, 
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code", 
     "collapsed": false, 
     "input": [
      "#!/usr/bin/python", 
      "\"\"\"", 
      "Demostration of the Hardy-Cross relaxation solution of the pipeline network", 
      " (1)  (2)   (3)    (7)", 
      "----o-----o-----o------", 
      "     \\(4) |(5) /(6)", 
      "      \\   |   /", 
      "       \\--o--/", 
      "       ", 
      "       ", 
      "\"\"\"", 
      "import numpy as n", 
      "import pylab as p", 
      "", 
      "L = n.array([2,2,3,3,2,3,2],dtype='f') # km", 
      "D = n.array([25,25,20,20,15,20,25],dtype='f') # cm", 
      "C = n.array([120,100,110,120,130,120,130],dtype='f') # Hazen Williams ", 
      "", 
      "# resistance per unit length as a function", 
      "r = lambda l,d,c: 1.526e7/(c**1.852*d**4.87)*l", 
      "# head loss as a function:", 
      "hf = lambda R,Q: R*n.sign(Q)*n.abs(Q)**1.852", 
      "", 
      "R = r(L,D,C)# resistances", 
      "", 
      "# define branches - each row contains numbers of pipes:", 
      "branch = n.array([[2,5,4],[3,6,5]])-1 # Python counts from zero, not 1", 
      "", 
      "rows,cols = branch.shape # rows = num of branches, cols = pipes in each", 
      "", 
      "", 
      "# initial guess that ", 
      "Q = n.array([500,250,250,-250,0,-250,500],dtype='f') # m**3/hr", 
      "", 
      "dQ = 500.0", 
      "", 
      "# main loop", 
      "while abs(dQ) > 0.1:", 
      "\tfor i in n.arange(rows):", 
      "\t\ty = hf(R,Q) ", 
      "\t\t", 
      "\t\tyq = n.abs(1.852*y/n.abs(Q))", 
      "\t\tyq[n.isnan(yq)] = 0.0", 
      "\t\t", 
      "\t\t# for the first branch", 
      "\t\t", 
      "\t\tsumyq = n.sum(yq[branch[i,:]])", 
      "\t\tsumy = n.sum(y[branch[i,:]])", 
      "\t\t", 
      "\t\tdQ = -1*sumy/sumyq", 
      "\t\tprint(\"dQ = %f\" % dQ)", 
      "\t\t", 
      "\t\tQ[branch[i,:]] += dQ", 
      "\t\t", 
      "", 
      "print(\"Discharges [m^3/hr]:\\n\")", 
      "print Q", 
      "", 
      "# we're looking for equivalent pipe solution with:", 
      "Deq = 25 # cm", 
      "Ceq = 120 # ", 
      "", 
      "y = hf(R[1],Q[1])+hf(R[2],Q[2])", 
      "", 
      "Leq = y/(r(1,25,120)*Q[0]**1.852)", 
      "", 
      "Leq = Leq + L[0] + L[6]*(120/C[6])**1.852", 
      "", 
      "print(\"Equivalent length = %3.2f km\" % Leq)", 
      "", 
      "", 
      "", 
      "", 
      "", 
      ""
     ], 
     "language": "python", 
     "outputs": [
      {
       "output_type": "stream", 
       "stream": "stdout", 
       "text": [
        "dQ = 70.284711", 
        "dQ = -18.083459", 
        "dQ = -6.477078", 
        "dQ = 0.765177", 
        "dQ = -0.338562", 
        "dQ = 0.067967", 
        "Discharges [m^3/hr]:", 
        "", 
        "[ 500.          313.46908569  232.74967957 -186.53092957   46.21876144", 
        " -267.25033569  500.        ]", 
        "Equivalent length = 7.44 km"
       ]
      }
     ], 
     "prompt_number": 6
    }, 
    {
     "cell_type": "code", 
     "collapsed": true, 
     "input": [], 
     "language": "python", 
     "outputs": []
    }
   ]
  }
 ]
}