cpipe

cpipe.ipynb

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()

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

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