Hardy cross method of solution for the pipeline network problem

Hardy Cross Method-Changed Github.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Resistances:  [0.00093889 0.00349946 0.00297863 0.0069502  0.00297863]\n",
      "common branch= 3\n",
      "dQ = -77.877620\n",
      "dQ = -52.802153\n",
      "dQ = -5.188758\n",
      "dQ = -4.059188\n",
      "dQ = -0.443181\n",
      "dQ = -0.045024\n",
      "Discharges [m^3/hr]:\n",
      "\n",
      "[-333.50955934  166.49044066  -26.60319425  193.0936349  -306.9063651 ]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\ADREW\\adrew-py\\lib\\site-packages\\ipykernel_launcher.py:98: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Implementation of the Hardy-Cross relaxation solution of the pipeline network, \n",
    "given in the following example:\n",
    "     (2)\n",
    "A------------B\n",
    "|          / |\n",
    "|         /  |\n",
    "|        /   |\n",
    "|(1) (3)/ (4)|\n",
    "|      /     |\n",
    "|     /      |\n",
    "|    /       |\n",
    "|   /        |\n",
    "|  /         |\n",
    "| /    (5)   |\n",
    "C------------D\n",
    "\n",
    "QA = 500 # m**3/h\n",
    "QB = 0   # m**3/h\n",
    "QC = 0   # m**3/h\n",
    "QD = -500 # m**3/h\n",
    "\n",
    "R1 = 0.000298\n",
    "R2 = 0.000939\n",
    "R3 = 0.00695\n",
    "R4 = 0.000349\n",
    "R5 = 0.000298\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "# import all the numerical library, replaces Matlab\n",
    "from numpy import * \n",
    "\n",
    "\n",
    "# Initial conditions:\n",
    "QA = 500 # m**3/h\n",
    "QB = 0   # m**3/h\n",
    "QC = 0   # m**3/h\n",
    "QD = -500 # m**3/h\n",
    "\n",
    "# array is like a vector in Matlab\n",
    "L = array([2,3,3,2,3],dtype='f') # Length in km, 'f' means 'float'\n",
    "D = array([25,20,20,15,20],dtype='f') # Diameter, cm\n",
    "C = 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\n",
    "r = lambda L,D,C: L*1.526e7/(C**1.852*D**4.87)\n",
    "\n",
    "R = r(L,D,C)\n",
    "print (\"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",
    "\n",
    "hf = lambda R,Q: R*sign(Q)*power(abs(Q),1.852)\n",
    "\n",
    "\n",
    "\n",
    "# define branches - each row contains numbers of pipes:\n",
    "branch = array([[2,3,1],[3,4,5]])-1 # Python counts from zero, not 1, the first pipe will be (0)\n",
    "\n",
    "rows,cols = branch.shape # rows = num of branches, cols = pipes in each\n",
    "\n",
    "#find common branch\n",
    "i=0\n",
    "for j in arange(cols):\n",
    "    for k in arange(cols):\n",
    "        if branch[i,j]==branch[i+1,k]:\n",
    "            common=branch[i+1,k]\n",
    "            break        \n",
    "print('common branch=',common+1)#count from 0 so we add +1\n",
    "\n",
    "# initial guess that \n",
    "Q = array([-250, 250.0, 0.0, 250, -250.0]) # m**3/hr\n",
    "\n",
    "dQ = 1.0\n",
    "\n",
    "# main loop\n",
    "while 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",
    "        \n",
    "\t\tQ[branch[i,:]] += dQ# += is equal to Q = Q + dQ\n",
    "        \n",
    "\t\tQ[common] = -Q[common]#Q[2] = -Q[2]\n",
    "\n",
    "print(\"Discharges [m^3/hr]:\\n\")\n",
    "print (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)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}

hardy_cross_pipe_network.ipynb