3923.sagews

#!/usr/bin/python
# function cpipe(filename)
# The program:
#      program cpipe, should be named cpipe.m
# 
#     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, should be saved separately as rough.m
#
def rough(ed,rn):
	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        
#end
#end




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=1e-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

data = np.loadtxt(DATA+'input.dat',usecols = (0,))
	
p1,p2,v1,v2,z1,z2,ha,hr,hm,d,l,vis,e,sw,q = data

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

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 = %4f m \n', 6.8367346938775508) (' actual energy removed between points 1 and 2 = %4f m \n', 0.0) (' minor losses between points 1 and 2 = %4f m \n', 0.0) (' length of conduit = %4f m \n', 110.0) (' conduit material roughness = %4f mm \n', 0.045000000000000005) (' fluid viscosity = %4f m^2/s \n', 9.9999999999999995e-07) 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 --------------------------------------------------------

Traceback (most recent call last): # diameter. the international system of units is used. File "", line 1, in <module> File "/tmp/tmpFh2J5i/___code___.py", line 166, in <module> fid.write('flow rate will be {:4f} m^3/s \n'.format(q)) ValueError: zero length field name in format