Download (right-click, save target as ...) this page as a Jupyterlab notebook from: CHYD-1

CE 3305 Engineering Fluid Mechanics
Summer 2022 Computational Hydraulics Exercise Set 1

LAST NAME, FIRST NAME

R00000000

Purpose :

Apply computational thinking (ENGR 1330) principles and hydraulics/fluid mechanics principles to analyze flow distribution in a pipe network.

Assessment Criteria :

Completion, plausible solutions, use JupyterLab as a calculator.

Problem 1.

Figure 1 is a gravity-flow pipe network with water supplied from a fixed-grade reservoir (pool elevation 100 meters) connected to node N2. All pipes are ductile iron.

Figure 1. Network Layout (Plan View)

The pipe dimensions and node demands are shown in the tables below.

Pipe ID Length(m) Diameter(mm)
1 1,220 254
2 1,829 254
3 1,829 305
4 1,982 610
5 2,134 254
6 915 457
7 1,524 254
8 91 305
Node ID Elevation(m) Demand(liters/sec)
N1 51.8 31.5
N2 54.9 31.5
N3 50.3 31.5
N4 47.3 94.6
N5 45.7 63.1
N6 44.2 94.6

Determine:

  1. The flow rate (and direction of flow) for each pipe in the network, for the case where the total head at the supply reservoir is 100 meters.
  2. The resultant pressure in SI units at each node.
  3. The Darcy-Weisbach friction factor for each ductile iron pipe of the network.
  4. The head loss from Node 2 to Node 6.
  5. The node with the lowest pressure.
In [ ]:
# sketch here
In [ ]:
# list known quantities
In [ ]:
# list unknown quantities

Input Data File

6
8
51.8 54.9 50.3 47.3 45.7 44.2
0.254 0.254 0.305 0.610 0.254 0.457 0.254 0.305
1220.0 1829.0 1829.0 1982.0 2134.0 915.0 1524.0 91.0
0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001
0.000011
1 1 1 1 1 1 1 1 1 1 1 1
1 0 -1 0 0 0 0 0
-1 -1 0 -1 0 0 0 1
0 1 0 0 -1 0 0 0
0 0 1 0 0 0 -1 0
0 0 0 1 0 1 1 0
0 0 0 0 1 -1 0 0
0.0315 0.0315 0.0315 0.0946 0.0631 0.0946 0 0 0 0 0 0 0 -100

Modify script for proper output labeling assuming SI inputs

In [ ]:
# governing principles (fluid mechanics)
In [ ]:
# solution algorithm
In [15]:
import numpy

def writeM(M,ir,jc,label):
    print ("------",label,"------")
    for i in range(0,ir,1):
        print (M[i][0:jc])
    print ("-----------------------------")
    return()

def writeV(V,ir,label):
    print ("------",label,"------")
    for i in range(0,ir,1):
        print (V[i])
    print ("-----------------------------")
    return()

def matrixmatrixmult(amatrix,bmatrix,rowNumA,colNumA,rowNumB,colNumB):
    AB =[[0.0 for j in range(colNumB)] for i in range(rowNumA)]
    for i in range(0,rowNumA):
        for j in range(0,colNumB):
            for k in range(0,colNumA):
                AB[i][j]=AB[i][j]+amatrix[i][k]*bmatrix[k][j]
    return(AB)

def matrixvectormult(amatrix,xvector,rowNumA,colNumA):
    bvector=[0.0 for i in range(rowNumA)]
    for i in range(0,rowNumA):
        for j in range(0,1):
            for k in range(0,colNumA):
                bvector[i]=bvector[i]+amatrix[i][k]*xvector[k]
    return(bvector)

def vectoradd(avector,bvector,length):
    cvector=[]
    for i in range(length):
        cvector.append(avector[i]+bvector[i])
    return(cvector)

def vectorsub(avector,bvector,length):
    cvector=[]
    for i in range(length):
        cvector.append(avector[i]-bvector[i])
    return(cvector)
             
def vdotv(avector,bvector,length):
    adotb=0.0
    for i in range(length):
        adotb=adotb+avector[i]*bvector[i]
    return(adotb)
# hydraulic elements prototype functions
# Jain Friction Factor Function -- Tested OK 23SEP16
import math   # This will import math module

def friction_factor(roughness,diameter,reynolds):
    temp1 = roughness/(3.7*diameter)
    temp2 = 5.74/(reynolds**(0.9))
    temp3 = math.log10(temp1+temp2)
    temp3 = temp3**2
    friction_factor = 0.25/temp3
    return(friction_factor)

# Velocity Function
def velocity(diameter,discharge):
    velocity=discharge/(0.25*math.pi*diameter**2)
    return(velocity)

# Reynolds Number Function  
def reynolds_number(velocity,diameter,mu):
    reynolds_number = abs(velocity)*diameter/mu
    return(reynolds_number)

# Geometric factor function
def k_factor(howlong,diameter,gravity):
    k_factor = (16*howlong)/(2.0*gravity*math.pi**2*diameter**5)
    return(k_factor)
# solution (using JupyterLab notebook) (computational thinking/algorithm development)
########## Initial Memory Allocations ###############
bvector = []
rowNumA = 0
colNumA = 0
rowNumB = 0
verbose = 'false' # set to true for in-class demonstration
#############################################

elevation = [] # null list node elevations
diameter =  [] # null list pipe diameters
distance =  [] # null list pipe lengths
roughness = [] # null list pipe roughness
flowguess = [] # null list pipe flow rates
nodearcs =  [] # node-arc incidence matrix
rhs_true =  [] # null list for nodal demands
tempvect = [] # null list for reading from external file, then recasting into one of the above lists

##############################################
# connect and read file for Pipeline Network #
##############################################
#infilename="L15-E2.txt"
infilename="chyd-1.txt"
afile = open(infilename,"r")
#afile = open("PipeNetwork.txt","r")
nnodes = int(afile.readline())
npipes = int(afile.readline())
# read elevation vector
tempvect.append([float(n) for n in afile.readline().strip().split()])
for i in range(0,nnodes,1):
    elevation.append(float(tempvect[0][i]))
tempvect = [] # reset vector
# read diameter vector
tempvect.append([float(n) for n in afile.readline().strip().split()])
for i in range(0,npipes,1):
    diameter.append(float(tempvect[0][i]))
tempvect = [] # reset vector
# read length vector
tempvect.append([float(n) for n in afile.readline().strip().split()])
for i in range(0,npipes,1):
    distance.append(float(tempvect[0][i]))
tempvect = [] # reset vector
# read roughness vector
tempvect.append([float(n) for n in afile.readline().strip().split()])
for i in range(0,npipes,1):
    roughness.append(float(tempvect[0][i]))
tempvect = [] # reset vector
# read viscosity (scalar)
viscosity = float(afile.readline())
# read current flow guess
tempvect.append([float(n) for n in afile.readline().strip().split()])
for i in range(0,npipes,1):
    flowguess.append(float(tempvect[0][i]))
tempvect = [] # reset vector
# read nodearc incidence matrix
## future revisions read directly into augmented matrix, or find way to release nodearc from stack
for irow in range(0,nnodes,1):  # then read each row
    nodearcs.append([float(n) for n in afile.readline().strip().split()])
# read demands guess
tempvect.append([float(n) for n in afile.readline().strip().split()])
for i in range(0,nnodes+npipes,1):
    rhs_true.append(float(tempvect[0][i]))
tempvect = [] # reset vector      
######################################
# end file read ,disconnect file     #
######################################
afile.close() # Disconnect the file

######################################
# echo the input in human readable   #
######################################
print('####ECHO INPUT################\nInput File: ',infilename)
print('number of nodes : ',nnodes)
print('number of pipes : ',npipes)
print('viscosity       : ',viscosity)
print ("-----------------------------")
for irow in range(0,nnodes):
    print('node id:',irow, ', elevation :',elevation[irow],' head :',rhs_true[irow+npipes])
print ("-----------------------------")
for jcol in range(0,npipes):
    print('pipe id:',jcol,', diameter : ' ,diameter[jcol],', distance : ',distance[jcol],
          ', roughness : ',roughness[jcol],', flow  : ',flowguess[jcol])
print ("-----------------------------")
##for jcol in range(0,nnodes+npipes):
##    print('irow :',jcol,' RHS True :',rhs_true[jcol])
##print ("-----------------------------")
print("node-arc incidence matrix")
for i in range(0,nnodes,1):
    print (nodearcs[i][0:npipes])
print ("-----------------------------")

# create augmented matrix
colNumA = npipes+nnodes
rowNumA = nnodes+npipes
augmentedMat = [] # null list to store augmented matrix

#######################################################################################
augmentedMat = [[0.0 for j in range(colNumA)]for i in range(rowNumA)] #fill with zeroes
#build upper left partition -- from nodearcs
for ir in range(0,nnodes):
    for jc in range (0,npipes):
        augmentedMat[ir][jc] = nodearcs[ir][jc]
istart=nnodes
iend=nnodes+npipes
jstart=npipes
jend=npipes+nnodes
for ir in range(istart,iend):
    for jc in range (jstart,jend):
        augmentedMat[ir][jc] = -1.0*nodearcs[jc-jstart][ir-istart] + 0.0
if verbose == 'true' :
    print("augmented matrix before loss factors")
    writeM(augmentedMat,rowNumA,colNumA,"augmented matrix")
    
#########Simulation Constants and Additional Memory Allocation #######
howmany=100 #iterations max
tolerance1 = 1e-27
tolerance2 = 1e-24
velocity_pipe = [0 for i in range(npipes)]  # null list velocities
reynolds      = [0 for i in range(npipes)]  # null list reynolds numbers
friction      = [0 for i in range(npipes)]  # null list friction 
geometry      = [0 for i in range(npipes)]  # null list geometry
lossfactor    = [0 for i in range(npipes)]  # null list loss
jacbMat = [] # null list to store jacobian matrix
jacbMat = [[0.0 for j in range(colNumA)]for i in range(rowNumA)] #fill with zeroes


solvecguess =[ 0.0 for i in range(rowNumA)] 
solvecnew =[ 0.0 for i in range(rowNumA)]
for i in range(0,npipes,1):
    solvecguess[i] = flowguess[i]
    geometry[i] = k_factor(distance[i],diameter[i],32.2)
#solvecguess is a current guess -- wonder if more pythonic way for this assignment
##    print('irow :',i,' Geometry Factor :',geometry[i])
##print ("-----------------------------")

###############################################################
## ITERATION LOOP                                             #
###############################################################
for iteration in range(howmany): # iteration outer loop
    if verbose == 'true' :
        print("solutions at begin of iteration",iteration)
        for jcol in range(0,nnodes+npipes):
            print('irow :',jcol,' solvecnew :',solvecnew[jcol]," solvecguess ",solvecguess[jcol])
        print ("-----------------------------")
    for i in range(0,npipes,1):
        velocity_pipe[i] = velocity(diameter[i],flowguess[i])    
        reynolds[i]=reynolds_number(velocity_pipe[i],diameter[i],viscosity)
        friction[i]=friction_factor(roughness[i],diameter[i],reynolds[i])
        lossfactor[i]=friction[i]*geometry[i]*abs(flowguess[i])
    if verbose == 'true' :
        for jcol in range(0,npipes):
            print('pipe id:',jcol,', velocity : ' ,velocity_pipe[jcol],', reynolds : ',reynolds[jcol],
          ', friction : ',friction[jcol],', loss factor : ',lossfactor[jcol],'flow guess',flowguess[jcol])
################################################################
# BUILD AUGMENTED MATRIX CURRENT Q+H SOLUTION                  #
################################################################
    augmentedMat = [[0.0 for j in range(colNumA)]for i in range(rowNumA)] #fill with zeroes
    #build upper left partition -- from nodearcs
    for ir in range(0,nnodes):
        for jc in range (0,npipes):
            augmentedMat[ir][jc] = nodearcs[ir][jc]
    #build lower right == transpose of upper left
    istart=nnodes
    iend=nnodes+npipes
    jstart=npipes
    jend=npipes+nnodes
    for ir in range(istart,iend):
        for jc in range (jstart,jend):
            augmentedMat[ir][jc] = -1.0*nodearcs[jc-jstart][ir-istart] + 0.0
    # build lower left partition of the matrix
    istart = nnodes
    iend = nnodes+npipes
    jstart = 0
    jend = npipes
    for i in range(istart,iend ):
        for j in range(jstart,jend ):
    #        print('i =',i,'j=',j)
            if (i-istart) == j :
    #            print('i =',i,'j=',j)
                augmentedMat[i][j] = -1.0*lossfactor[j] + 0.0
    if verbose == 'true' :
        print("updated augmented matrix in iteration",iteration)
        writeM(augmentedMat,rowNumA,colNumA,"augmented matrix")
################################################################
# BUILD JACOBIAN MATRIX CURRENT Q+H SOLUTION                   #
################################################################        
    # now build current jacobian
    for i in range(rowNumA):
        for j in range(colNumA):
            jacbMat[i][j] = augmentedMat[i][j]
    # modify lower left partition
    istart = nnodes
    iend = nnodes+npipes
    jstart = 0
    jend = npipes
    for i in range(istart,iend ):
        for j in range(jstart,jend ):
    #        print('i =',i,'j=',j)
            if (i-istart) == j :
    #            print('i =',i,'j=',j)
                jacbMat[i][j] = 2.0*jacbMat[i][j]
##        for jcol in range(0,nnodes+npipes):
##            print('irow :',jcol,' solvecnew :',solvecnew[jcol]," solvecguess ",solvecguess[jcol])
##        print ("-----------------------------")
# matrix multiply augmentedMat*solvecguess to get current g(Q)
#    gq = [0.0 for i in range(rowNumA)] # zero gradient vector
##    if verbose == 'true' :
##        print("augmented matrix in iteration",iteration)
##        writeM(augmentedMat,rowNumA,colNumA,"augmented matrix before mmult")
    gq = matrixvectormult(augmentedMat,solvecguess,rowNumA,colNumA)
##    if verbose == 'true' :
##        writeV(gq,rowNumA,"gq vectorbefore subtract rhs_true")
# subtract rhs
#    for i in range(rowNumA):
    gq = vectorsub(gq,rhs_true,rowNumA)#vector subtract
    if verbose == 'true' :
        print("computed g(q) in iteration",iteration)
        writeV(gq,rowNumA,"gq vector")
        print("compare current and new guess")
        for jcol in range(0,nnodes+npipes):
            print('irow :',jcol,' solvecnew :',solvecnew[jcol]," solvecguess ",solvecguess[jcol])
        print ("-----------------------------")
    dq = [0.0 for i in range(rowNumA)] # zero update vector
    if verbose == 'true' :
        writeV(dq,rowNumA,"dq vector before linear solve")
    if verbose == 'true' :
        print("jacobian before linearsolve in iteration",iteration)
        writeM(jacbMat,rowNumA,colNumA,"jabobian matrix")
#    dq = linearsolver(jacbMat,gq) # memory leak after this call - linearsolve clobbers input lists
    dq = numpy.linalg.solve(jacbMat,gq) #the numpy equivalent
    if verbose == 'true' :
        print("jacobian after linearsolve in iteration",iteration)
        writeM(jacbMat,rowNumA,colNumA,"jabobian matrix")

    if verbose == 'true' :
        writeV(dq,rowNumA,"dq vector -after linear solve")
    solvecnew = vectorsub(solvecguess,dq,rowNumA)#vector subtract
    if verbose == 'true' :
        print("Q_new = Q_old - DQ")
        writeV(solvecnew,rowNumA,"new guess vector")
#    tempvect =[ 0.0 for i in range(rowNumA)]
##        tempvect = matrixvectormult(jacbMat,dq,rowNumA,colNumA)
##        writeV(tempvect,rowNumA,"J*dq vector")
##        tempvect = vectorsub(tempvect,gq,rowNumA)
##        writeV(tempvect,rowNumA,"J*dq - gq vector")
        print("just after computing new guess, should be different")
        for jcol in range(0,nnodes+npipes):
            print('irow :',jcol,' solvecnew :',solvecnew[jcol]," solvecguess ",solvecguess[jcol])
        print ("-----------------------------")
#test for stopping
    tempvect =[ 0.0 for i in range(rowNumA)]
    for i in range(rowNumA):
        tempvect[i] = abs(solvecnew[i] - solvecguess[i])
    test1 = vdotv(tempvect,tempvect,rowNumA)
    if verbose == 'true' :
        print('test1',test1)
    tempvect =[ 0.0 for i in range(rowNumA)]
    for i in range(rowNumA):
        tempvect[i] = abs(gq[i])
    test2 = vdotv(tempvect,tempvect,rowNumA)
    if verbose == 'true' :
        print('test2',test2)
    if test1 < tolerance1 :
        print("###EXIT TYPE###\nUpdate not changing --exit and report current update")
        print("iteration",iteration)
# update guess
        solvecguess[:] = solvecnew[:]
        for i in range(0,npipes,1):
            flowguess[i] = solvecguess[i]

        break
    if test2 < tolerance2 :
        print("###EXIT TYPE###\nGradient near zero --exit and report current update")
        print("iteration",iteration)
# update guess
        solvecguess[:] = solvecnew[:]
        for i in range(0,npipes,1):
            flowguess[i] = solvecguess[i]

        break
    if verbose == 'true' :
        print("solution continuing")
        print("iteration",iteration)
    # update guess
    solvecguess[:] = solvecnew[:]
    if verbose == 'true' :
        for i in range(0,npipes,1):
            flowguess[i] = solvecguess[i]
## Write Current State ######################
        gq = matrixvectormult(augmentedMat,solvecguess,rowNumA,colNumA)
        print('number of nodes : ',nnodes)
        print('number of pipes : ',npipes)
        print('viscosity       : ',viscosity)
        print ("-----------------------------")
        for irow in range(0,nnodes):
            print('node id:',irow, ', elevation :',elevation[irow])
        print ("-----------------------------")
        for jcol in range(0,npipes):
            print('pipe id:',jcol,', diameter : ' ,diameter[jcol],', distance : ',distance[jcol],
          ', roughness : ',roughness[jcol],', flow guess : ',round(flowguess[jcol],3))
        print ("-----------------------------")
        for jcol in range(0,nnodes+npipes):
            print('irow :',jcol,' RHS True :',rhs_true[jcol],"RHS Current",round(gq[jcol],3))
        print ("-----------------------------")
        for jcol in range(0,nnodes+npipes):
            print('irow :',jcol,' solvecnew :',solvecnew[jcol]," solvecguess ",solvecguess[jcol])
        print ("-----------------------------")   
################################################
# end of outer loop
################################################

# Report Results
print("#####SIMULATION RESULTS#####\nResults at iteration = :",iteration)
for i in range(0,npipes,1):
    flowguess[i] = solvecguess[i]
print('number of nodes : ',nnodes)
print('number of pipes : ',npipes)
print('viscosity       : ',viscosity)
print ("-----------------------------")
istart = int(npipes)
for irow in range(0,nnodes):
    print('node id:',irow+1, ', elevation (m) :',elevation[irow],' head (m) :',round(solvecnew[irow+npipes],3),' pressure (Pa) :', round((9800)*(solvecnew[irow+npipes]-elevation[irow]),3))
print ("-----------------------------")
for jcol in range(0,npipes):
    print('pipe id:',jcol+1,', diameter (m) : ' ,diameter[jcol],', distance (m) : ',distance[jcol],
  ', friction factor : ',round(friction[jcol],3),', flow (cms) : ',round(flowguess[jcol],3))
print ("-----------------------------")
if verbose == 'true' :
    for jcol in range(0,nnodes+npipes):
        print('irow :',jcol,' RHS True :',rhs_true[jcol],"RHS Current",gq[jcol])
    print ("-----------------------------")
    for jcol in range(0,nnodes+npipes):
        print('irow :',jcol,' solvecnew :',solvecnew[jcol]," solvecguess ",solvecguess[jcol])
    print ("-----------------------------")

####ECHO INPUT################
Input File:  chyd-1.txt
number of nodes :  6
number of pipes :  8
viscosity       :  1.1e-05
-----------------------------
node id: 0 , elevation : 51.8  head : 0.0
node id: 1 , elevation : 54.9  head : 0.0
node id: 2 , elevation : 50.3  head : 0.0
node id: 3 , elevation : 47.3  head : 0.0
node id: 4 , elevation : 45.7  head : 0.0
node id: 5 , elevation : 44.2  head : -100.0
-----------------------------
pipe id: 0 , diameter :  0.254 , distance :  1220.0 , roughness :  1e-05 , flow  :  1.0
pipe id: 1 , diameter :  0.254 , distance :  1829.0 , roughness :  1e-05 , flow  :  1.0
pipe id: 2 , diameter :  0.305 , distance :  1829.0 , roughness :  1e-05 , flow  :  1.0
pipe id: 3 , diameter :  0.61 , distance :  1982.0 , roughness :  1e-05 , flow  :  1.0
pipe id: 4 , diameter :  0.254 , distance :  2134.0 , roughness :  1e-05 , flow  :  1.0
pipe id: 5 , diameter :  0.457 , distance :  915.0 , roughness :  1e-05 , flow  :  1.0
pipe id: 6 , diameter :  0.254 , distance :  1524.0 , roughness :  1e-05 , flow  :  1.0
pipe id: 7 , diameter :  0.305 , distance :  91.0 , roughness :  1e-05 , flow  :  1.0
-----------------------------
node-arc incidence matrix
[1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[-1.0, -1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0]
[0.0, 1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0, 0.0]
[0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0]
-----------------------------
###EXIT TYPE###
Update not changing --exit and report current update
iteration 46
#####SIMULATION RESULTS#####
Results at iteration = : 46
number of nodes :  6
number of pipes :  8
viscosity       :  1.1e-05
-----------------------------
node id: 1 , elevation (m) : 51.8  head (m) : 70.068  pressure (Pa) : 179023.337
node id: 2 , elevation (m) : 54.9  head (m) : 95.728  pressure (Pa) : 400118.228
node id: 3 , elevation (m) : 50.3  head (m) : 83.623  pressure (Pa) : 326570.271
node id: 4 , elevation (m) : 47.3  head (m) : 62.109  pressure (Pa) : 145128.174
node id: 5 , elevation (m) : 45.7  head (m) : 93.561  pressure (Pa) : 469040.665
node id: 6 , elevation (m) : 44.2  head (m) : 91.714  pressure (Pa) : 465637.732
-----------------------------
pipe id: 1 , diameter (m) :  0.254 , distance (m) :  1220.0 , friction factor :  0.014 , flow (cms) :  0.064
pipe id: 2 , diameter (m) :  0.254 , distance (m) :  1829.0 , friction factor :  0.014 , flow (cms) :  0.02
pipe id: 3 , diameter (m) :  0.305 , distance (m) :  1829.0 , friction factor :  0.014 , flow (cms) :  0.032
pipe id: 4 , diameter (m) :  0.61 , distance (m) :  1982.0 , friction factor :  0.016 , flow (cms) :  0.232
pipe id: 5 , diameter (m) :  0.254 , distance (m) :  2134.0 , friction factor :  0.014 , flow (cms) :  -0.011
pipe id: 6 , diameter (m) :  0.457 , distance (m) :  915.0 , friction factor :  0.015 , flow (cms) :  -0.106
pipe id: 7 , diameter (m) :  0.254 , distance (m) :  1524.0 , friction factor :  0.014 , flow (cms) :  -0.062
pipe id: 8 , diameter (m) :  0.305 , distance (m) :  91.0 , friction factor :  0.014 , flow (cms) :  0.347
-----------------------------
In [17]:
# discussion
In [ ]: