Hardy cross method of solution for the pipeline network problem

hardy_cross_pipe_network.ipynb

{
 "metadata": {
  "name": "hardy_cross_pipe_network"
 }, 
 "name": "hardy_cross_pipe_network", 
 "nbformat": 2, 
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code", 
     "collapsed": false, 
     "input": "\"\"\"\nImplementation of the Hardy-Cross relaxation solution of the pipeline network, \ngiven in the following example:\n     (2)\nA------------B\n|          / |\n|         /  |\n|        /   |\n|(1) (3)/ (4)|\n|      /     |\n|     /      |\n|    /       |\n|   /        |\n|  /         |\n| /    (5)   |\nC------------D\n\nQA = 500 # m**3/h\nQB = 0   # m**3/h\nQC = 0   # m**3/h\nQD = -500 # m**3/h\n\nR1 = 0.000298\nR2 = 0.000939\nR3 = 0.00695\nR4 = 0.000349\nR5 = 0.000298\n\n\"\"\"\n\n# import all the numerical library, replaces Matlab\nfrom numpy import * \n\n\n# Initial conditions:\nQA = 500 # m**3/h\nQB = 0   # m**3/h\nQC = 0   # m**3/h\nQD = -500 # m**3/h\n\n# array is like a vector in Matlab\nL = array([2,3,3,2,3],dtype='f') # Length in km, 'f' means 'float'\nD = array([25,20,20,15,20],dtype='f') # Diameter, cm\nC = array([100,110,120,130,120],dtype='f') # Hazen Williams friction coefficient\n \n\n# resistance as a single line function that receives C, D, L and returns R\nr = lambda L,D,C: L*1.526e7/(C**1.852*D**4.87)\n\nR = r(L,D,C)\nprint \"Resistances: \", R\n\n# Multiline definition of a function\n# starts with def name of the function and inputs in parentheses\n# next line is with indentation: 4 spaces or a TAB\n# ends with the \"return\" and the list of outputs\n# see the following example\n\n#def resistance(L,d,c):\n#    r = 1.526e7/(c**1.852*d**4.87)*L\n#    return r           \n\n\n# head loss as a function\n# note that one cannot raise a negative value to the power of 1.852\n\nhf = lambda R,Q: R*sign(Q)*power(abs(Q),1.852)\n\n\n\n# define branches - each row contains numbers of pipes:\nbranch = array([[2,3,1],[3,4,5]])-1 # Python counts from zero, not 1, the first pipe will be (0)\n\nrows,cols = branch.shape # rows = num of branches, cols = pipes in each\n\n\n# initial guess that \nQ = array([-250, 250.0, 0.0, 250, -250.0]) # m**3/hr\n\ndQ = 1.0\n\n# main loop\nwhile abs(dQ) > 0.5:\n    # arange(N) gives a vector from 0 to N-1, i.e. arange(3) = 0,1,2\n\tfor i in arange(rows):\n                # estimate the losses in each pipe\n\t\ty = hf(R,Q) \n\t\t\n\n\t\tyq = abs(1.852*y/abs(Q))\n                \n                # Remove NaNs for the exceptional cases\n\t\tyq[isnan(yq)] = 0.0\n\t\t\n\t\t# for the first branch\n\t\t# Sum is a sum :)\n\t\tsumyq = sum(yq[branch[i,:]])\n\t\tsumy = sum(y[branch[i,:]])\n\t\t\n                # Estimate the correction dQ for the next step:\n\t\tdQ = -1*sumy/sumyq\n\n\t\tprint(\"dQ = %f\" % dQ)\n\t\t\n                # += is equal to Q = Q + dQ\n\t\tQ[branch[i,:]] += dQ\n\t\t\n\nprint(\"Discharges [m^3/hr]:\\n\")\nprint Q\n\n# we're looking for equivalent pipe solution with:\n# Deq = 25 # cm\n# Ceq = 120 # \n\n# y = hf(R[1],Q[1])+hf(R[2],Q[2])\n\n# Leq = y/(r(1,25,120)*Q[0]**1.852)\n\n# Leq = Leq + L[0] + L[6]*(120/C[6])**1.852\n\n# print(\"Equivalent length = %3.2f km\" % Leq)", 
     "language": "python", 
     "outputs": [
      {
       "output_type": "stream", 
       "stream": "stdout", 
       "text": "Resistances:  [ 0.00093889  0.00349946  0.00297863  0.0069502   0.00297863]\ndQ = -77.877620\ndQ = -44.395536\ndQ = 14.902768\ndQ = -2.315951\ndQ = 0.672938\ndQ = -0.092987\nDischarges [m^3/hr]:\n\n[-312.30191449  187.69808551 -109.10638868  203.19552581 -296.80447419]"
      }, 
      {
       "output_type": "stream", 
       "stream": "stderr", 
       "text": "-c:90: RuntimeWarning: invalid value encountered in divide\n-c:90: RuntimeWarning: invalid value encountered in absolute"
      }
     ], 
     "prompt_number": 1
    }
   ]
  }
 ]
}