Skip to content

Instantly share code, notes, and snippets.

@alexlib

alexlib/cpipe.ipynb

Created November 17, 2012 20:38
Show Gist options
  • Save alexlib/4100013 to your computer and use it in GitHub Desktop.
Save alexlib/4100013 to your computer and use it in GitHub Desktop.
Download ZIP
cpipe
Display the source blob
Display the rendered blob
Raw
cpipe.ipynb
{
"metadata": {
"name": "cpipe"
},
"nbformat": 2,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"source": [
"## Pipe flow - estimate of pressure drop due to friction",
"### For the following input data variations, please estimate the discharge. ",
"",
"*input data:*",
"",
" p1\t- gauge pressure at point 1 in kPa:\t\t0.0 ",
" p2\t- gauge pressure at point 2 in kPa:\t\t0.0 ",
" v1\t- velocity if essentially zero, otherwise 1.0:\t0.0 ",
" v2\t- velocity if essentially zero, otherwise 1.0:\t1.0 ",
" z1\t- elevation at point 1 in meters:\t\t\t80.0 ",
" z2\t- elevation at point 2 in meters:\t\t\t60.; 50.; 40.; 30.; 20.; ",
" ha\t- energy added between points 1 and 2 in kilowatts:\t0.; 50.; 100.;\t ",
" hr\t- energy removed between points 1 and 2 in kilowatts:\t0. ",
" hm\t- minor losses between points 1 and 2 in meters:\t0.; ",
" d\t- diameter of a pipe in millimeters: \t\t\tfrom 300. \u2013 to 400. (delta=20 mm) ",
" l\t- length of a pipe in meters: \t\t\t\t1000 ",
" vis\t- kinematic viscosity in m^2/s:\t\t\t0.406e-6 ",
" e\t- pipe material roughness in meters; zero for smooth: 0.0005 ",
" sw\t- specific weight of a fluid in kilonewtons/m^3:\t7.05 ",
" q\t- flow rate in m^3/s:\t\t\t\t\t0. (unknown!) "
]
},
{
"cell_type": "markdown",
"source": [
"## Exercise no. 3:",
"1. For each $h_a$ (see above), for each $z_1 - z_2$ (see above), and for each diameter in the given range of $d$, estimate the discharge $q$",
"2. The result should be prepared as a figure of $q(d)$ ",
"",
"## The code:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"'''",
"function cpipe()",
"",
"The program:",
" program cpipe, renamed from cpipe.m",
" ",
" The function is not written very well, resembles the older Fortran/Matlab version",
"",
" this program computes the flow rate or the required conduit",
" diameter. the international system of units is used. ",
" the application of the program is limited to cases involving ",
" a single pipe with constant diameter.",
" ",
" this program is based on solution of the energy equation. ",
" accordingly, certain data must be entered for each of two points ",
" in a pipe system. however, both the velocity at point 1 and the ",
" velocity at point 2 are not considered as known values initially. ",
" if either value of velocity is essentially zero (such as at a ",
" reservoir or tank surface) - enter 0 (zero) for that velocity or ",
" for both velocities if both are essentially zero, otherwise - enter",
" 1.0 for all other velocities. ",
"'''",
"",
"import numpy as np",
"from numpy import pi,sqrt,log10",
"import sys",
"import os",
"",
"",
"#####################################################################################",
"# ",
"# function \"rough\"",
"#",
"def rough(ed,rn):",
"\t''' rough(ed,rn) estimates the friction coefficient f",
"\tInputs: ed = epsilon/D = relative roughness, dimensionless",
"\t\t rn = Reynolds number",
"\tOutput: f = friction coefficient according to Colebrook-White equation",
"\t'''",
"\tif rn < 2000: # if it's a laminar flow",
"\t\tf = 64/rn",
"\t\treturn f",
"\t",
" # friction coefficient quess",
"\tf = 0.006",
"\ttry1 = 1./np.sqrt(f)+2*np.log10(ed/3.7+2.51/rn/np.sqrt(f))",
"\tfor i in np.arange(10000):",
"\t\tf = f + 0.00001",
"\t\ttry2 = 1./sqrt(f)+2.*log10(ed/3.7+2.51/rn/sqrt(f))",
"\t\tif try1*try2>0:",
"\t\t\ttry1 = try2",
"\t\telse:",
"\t\t\tf = f - 0.000005",
"",
"\treturn f ",
"",
"",
"factor = 1000 # factor to convert mm into m",
"#--------------------------------------------------------------------------",
"# Input Data",
"#",
"# ",
"",
"",
"# Try to read from the file, otherwise, use defaults:",
"try:",
"\tdata = np.loadtxt('input.dat',usecols = (0,))",
"\tp1,p2,v1,v2,z1,z2,ha,hr,hm,d,l,vis,e,sw,q = data",
"except:",
"\t# Defaults:",
"",
" p1 = 1.0 # gauge pressure at point 1 in kilopascals",
" p2 = 0.0 # gauge pressure at point 2 in kilopascals",
" v1 = 0.0 # velocity = if essentially zero, otherwise 1.0",
" v2 = 1.0 # velocity = if essentially zero, otherwise 1.0",
" z1 = 39.0 # elevation at point 1 in meters",
" z2 = 0.0 # elevation at point 2 in meters",
" ha = 0.0 # energy added between points 1 and 2 in kilowats",
" hr = 0.0 # energy removed between points 1 and 2 in kilowats",
" hm = 2.0 # minor losses between points 1 and 2 in meters",
" d = 100.0 # diameter of conduit in millimeters",
" l = 100.0 # length of conduit in meters",
" vis = 1.0e-6 # kinematic viscosity in m^2/s",
" e = 0.000045 # conduit material roughness in meters zero for smooth",
" sw = 9.8 # specific weight of fluid in kilonewtons per m^3",
" q = 0.0 # flow rate in m^3/s, =0 if unknown",
"\t",
"#--------------------------------------------------------------------------",
"g = 9.807 # gravity",
"p1sw = p1/sw",
"p2sw = p2/sw",
"",
"# friction coefficient initial guess",
"ff = 0.025",
"",
"# for i=1,11 #do 333",
"# d = 200 + (i-1)*50",
"",
"# d is known, find q",
"if v1 > 0.0001:",
"\tv1 = 1.0/2.0/g",
"",
"if v2>0.0001:",
"\tv2 = 1.0/2.0/g",
"",
"for j in range(1,10000):",
"\thf = ff*l/d*factor/2/g #105 ",
"\that = ha/sw",
"\thrt = hr/sw",
" # flow rate q first guess",
"\tq = 0.001",
"\tvtry = (q/(pi*(d/factor)**2/4))**2",
"\ttry1 = p1sw+vtry*v1+z1+hat/q-hrt/q-(p2sw+vtry*v2+z2+hf*vtry)",
"\tfor n in range(1,10000):",
"\t\tq = q + 0.001",
"\t\tvtry = (q/(pi*(d/factor)**2/4))**2",
"\t\ttry2 = p1sw+vtry*v1+z1+hat/q-hrt/q-(p2sw+vtry*v2+z2+hf*vtry)",
"\t\tif try1*try2>0:",
"\t\t\ttry1 = try2",
"\t\telse:",
"\t\t\tq = q - 0.0005",
"\t\t\tbreak",
"\t \t#end",
"\t#end",
"\tv = q/(pi*(d/factor)**2/4)",
"\trn = d/factor*v/vis",
"\ted = e/d*factor",
" # compute friction coefficient",
"\tf=rough(ed,rn) # call 'rough' subroutine",
"\tdiff = abs(f-ff)",
"\tif diff<0.0001: # goto 104",
"\t\tbreak",
"\telse:",
"\t\tff = f",
"\t\tcontinue #goto 105",
"\t#end",
"# end",
"if v1 >= 0.0001: # 104 ",
"\tv1 = v",
"",
"if v2 >= 0.0001:",
"\tv2 = v",
"",
" ",
"# printing the input / output",
"' \\n'",
"print 'sample analysis for a closed conduit carrying fluid \\n'",
"print '--------------------------------------------------------\\n'",
"print ''",
"print' input data: \\n'",
"print' -----------\\n'",
"print' \\n'",
"print' pressure at point 1 = %4f kPa \\n' % p1",
"print' pressure at point 2 = %4f kPa \\n' % p2",
"print' elevation at point 1 = %4f m \\n' % z1",
"print' elevation at point 2 = %4f m \\n' % z2",
"print' actual energy added between points 1 and 2 = %4f m \\n' % hat",
"print' actual energy removed between points 1 and 2 = %4f m \\n' % hrt",
"print' minor losses between points 1 and 2 = %4f m \\n' % hm",
"print' length of conduit = %4f m \\n' % l ",
"print' conduit material roughness = %4f mm \\n' % (e*factor)",
"print' fluid viscosity = %4f m^2/s \\n' % vis",
"print' \\n'",
"print' computed results: \\n'",
"print' ---------------\\n'",
"print' \\n'",
"print' flow rate will be %4f m^3/s \\n' % q",
"print' velocity at point 1 = %4f m/s \\n' % v1",
"print' velocity at point 2 = %4f m/s \\n' % v2",
"print'--------------------------------------------------------\\n'",
"",
"# writing output in a file: out.txt",
"fid = open('out.txt','w')",
"fid.write('computed data: \\n')",
"fid.write('flow rate will be {:4f} m^3/s \\n'.format(q))",
"fid.write('velocity at point 1 = {:4f} m/s \\n'.format(v1))",
"fid.write('velocity at point 2 = {:4f} m/s \\n'.format(v2))",
"fid.close()",
""
],
"language": "python",
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"sample analysis for a closed conduit carrying fluid ",
"",
"--------------------------------------------------------",
"",
"",
" input data: ",
"",
" -----------",
"",
" ",
"",
" pressure at point 1 = 0.000000 kPa ",
"",
" pressure at point 2 = 0.000000 kPa ",
"",
" elevation at point 1 = 39.000000 m ",
"",
" elevation at point 2 = 0.000000 m ",
"",
" actual energy added between points 1 and 2 = 6.836735 m ",
"",
" actual energy removed between points 1 and 2 = 0.000000 m ",
"",
" minor losses between points 1 and 2 = 0.000000 m ",
"",
" length of conduit = 110.000000 m ",
"",
" conduit material roughness = 0.045000 mm ",
"",
" fluid viscosity = 0.000001 m^2/s ",
"",
" ",
"",
" computed results: ",
"",
" ---------------",
"",
" ",
"",
" flow rate will be 0.477500 m^3/s ",
"",
" velocity at point 1 = 0.000000 m/s ",
"",
" velocity at point 2 = 6.755243 m/s ",
"",
"--------------------------------------------------------"
]
}
],
"prompt_number": 5
},
{
"cell_type": "code",
"collapsed": true,
"input": [],
"language": "python",
"outputs": [],
"prompt_number": 5
}
]
}
]
}
Raw
cpipe_alex.py
#!/usr/bin/python
'''
function cpipe(filename)
The program:
program cpipe, renamed from cpipe.m
The function is not written very well, resembles the older Fortran/Matlab version
this program computes the flow rate or the required conduit
diameter. the international system of units is used.
the application of the program is limited to cases involving
a single pipe with constant diameter.
this program is based on solution of the bernoulli equation.
accordingly, certain data must be entered for each of two points
in a pipe system. however, both the velocity at point 1 and the
velocity at point 2 are not considered as known values initially.
if either value of velocity is essentially zero (such as at a
reservoir or tank surface) - enter 0 (zero) for that velocity or
for both velocities if both are essentially zero, otherwise - enter
1.0 for all other velocities.
'''
import numpy as np
from numpy import pi,sqrt,log10
import sys
import os
#####################################################################################
#
# function "rough"
#
def rough(ed,rn):
''' rough(ed,rn) estimates the friction coefficient f
Inputs: ed = epsilon/D = relative roughness, dimensionless
rn = Reynolds number
Output: f = friction coefficient according to Colebrook-White equation
'''
if rn<2000:
f = 64/rn
return f
#end
# friction coefficient quess
f = 0.006
try1 = 1./np.sqrt(f)+2*np.log10(ed/3.7+2.51/rn/np.sqrt(f))
for i in np.arange(10000):
f = f + 0.00001
try2 = 1./sqrt(f)+2.*log10(ed/3.7+2.51/rn/sqrt(f))
if try1*try2>0:
try1 = try2
else:
f = f - 0.000005
return f
factor = 1000 # factor to convert mm into m
#--------------------------------------------------------------------------
# Input Data
#
#
# Defaults:
p1 = 0 # gauge pressure at point 1 in kilopascals
p2 = 0 # gauge pressure at point 2 in kilopascals
v1 = 0 # velocity = if essentially zero, otherwise 1.0
v2 = 1 # velocity = if essentially zero, otherwise 1.0
z1 =39 # elevation at point 1 in meters
z2 = 0 # elevation at point 2 in meters
ha = 0 # energy added between points 1 and 2 in kilowats
hr = 0 # energy removed between points 1 and 2 in kilowats
hm = 2 # minor losses between points 1 and 2 in meters
d = 100 # diameter of conduit in millimeters
l = 100 # length of conduit in meters
vis = 1.0e-6 # kinematic viscosity in m^2/s
e = 0.000045 # conduit material roughness in meters zero for smooth
sw = 9.8 # specific weight of fluid in kilonewtons per m^3
q = 0 # flow rate in m^3/s, =0 if unknown
# Try to read from the file, otherwise, use defaults:
try:
data = np.loadtxt('input.dat',usecols = (0,))
p1,p2,v1,v2,z1,z2,ha,hr,hm,d,l,vis,e,sw,q = data
except:
pass
#--------------------------------------------------------------------------
g = 9.807
p1sw = p1/sw
p2sw = p2/sw
# friction coefficient initial guess
ff = 0.02
# for i=1,11 #do 333
# d = 200 + (i-1)*50
# d is known, find q
if v1>0.0001:
v1 = 1/2/g
if v2>0.0001:
v2 = 1./2./g
for j in range(1,10000):
hf = ff*l/d*factor/2/g #105
hat = ha/sw
hrt = hr/sw
# flow rate q first guess
q = 0.001
vtry = (q/(pi*(d/factor)**2/4))**2
try1 = p1sw+vtry*v1+z1+hat/q-hrt/q-(p2sw+vtry*v2+z2+hf*vtry)
for n in range(1,10000):
q = q + 0.001
vtry = (q/(pi*(d/factor)**2/4))**2
try2 = p1sw+vtry*v1+z1+hat/q-hrt/q-(p2sw+vtry*v2+z2+hf*vtry)
if try1*try2>0:
try1 = try2
else:
q = q - 0.0005
break
#end
#end
v = q/(pi*(d/factor)**2/4)
rn = d/factor*v/vis
ed = e/d*factor
# compute friction coefficient
f=rough(ed,rn) # call rough
diff = abs(f-ff)
if diff<0.0001: # goto 104
break
else:
ff = f
continue #goto 105
#end
# end
if v1>=0.0001: # 104
v1 = v
#end
if v2>=0.0001:
v2 = v
# end
# writing input and output on screen
(' \n')
print(' sample analysis for a closed conduit carrying fluid \n')
print('--------------------------------------------------------\n')
print(' \n')
print(' input data: \n')
print(' \n')
print(' pressure at point 1 = %4f kPa \n' % p1)
print(' pressure at point 2 = %4f kPa \n' % p2)
print(' elevation at point 1 = %4f m \n' % z1)
print(' elevation at point 2 = %4f m \n' % z2)
print(' actual energy added between points 1 and 2 = %4f m \n',hat)
print(' actual energy removed between points 1 and 2 = %4f m \n',hrt)
print(' minor losses between points 1 and 2 = %4f m \n',hm)
print(' length of conduit = %4f m \n',l)
print(' conduit material roughness = %4f mm \n',e*factor)
print(' fluid viscosity = %4f m^2/s \n',vis)
print(' \n')
print(' computed data: \n')
print(' \n')
print(' flow rate will be %4f m^3/s \n' % q)
print(' velocity at point 1 = %4f m/s \n' % v1)
print(' velocity at point 2 = %4f m/s \n' % v2)
print('--------------------------------------------------------\n')
# writing output in a file: out.txt
fid = open('out.txt','w')
fid.write('computed data: \n')
fid.write('flow rate will be {:4f} m^3/s \n'.format(q))
fid.write('velocity at point 1 = {:4f} m/s \n'.format(v1))
fid.write('velocity at point 2 = {:4f} m/s \n'.format(v2))
fid.close()
Raw
input.dat
0.0 p1 - gauge pressure at point 1 in kilopascals
0.0 p2 - gauge pressure at point 2 in kilopascals
0.0 v1 - velocity - 0 if essentially zero, otherwise 1.0
1.0 v2 - velocity - 0 if essentially zero, otherwise 1.0
39.0 z1 - elevation at point 1 in meters
0.0 z2 - elevation at point 2 in meters
67.0 ha - energy added between points 1 and 2 in kilowats
0.0 hr - energy removed between points 1 and 2 in kilowats
0.0 hm - minor losses between points 1 and 2 in meters
300. d - diameter of conduit in millimeters
110.0 l - length of conduit in meters
1.0e-6 vis - kinematic viscosity in m^2/s
0.000045 e - conduit material roughness in meters; zero for smooth
9.8 sw - specific weight of fluid in kilonewtons per m^3
0.0 q - flow rate in m^3/s, =0 if unknown
Raw
out.txt
computed data:
flow rate will be 0.477500 m^3/s
velocity at point 1 = 0.000000 m/s
velocity at point 2 = 6.755243 m/s
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment