/
drain3.py
88 lines (66 loc) · 2.02 KB
/
drain3.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
# -*- coding: utf-8 -*-
"""
Time to drain the water from a horizontal pipe through an orifice at the bottom.
Jet speed is v = sqrt(2gh), being h the height of water inside the pipe.
chord = diam * sqrt(1 - ((h - r)/r)**2)
Created on Tue May 19 15:46:53 2015
@author: rarossi
"""
from math import sin, asin, pi
from matplotlib import pyplot as plt
from scipy.optimize import bisect
g = 9.806 # m/s^2
_debug = False
def chord(r, d):
"""Return the length of the chord of a circle of radius r at a distance d
from the centre."""
return 2 * r * (1 - (d/r)**2)**0.5
def area_above_segment(r, h):
"""radius, height --> area of the circular segment above a chord at height"""
c = chord(r, abs(h-r))
theta = 2 * asin(c/2/r)
area_segment = r**2 / 2 * (theta - sin(theta))
if h >= r:
return area_segment
else:
return pi*r**2 - area_segment
def calc_dh(area, h0, r):
if _debug: print(area, h0)
def func(h):
return area_above_segment(r, h) - area_above_segment(r, h0) - area
return h0 - bisect(func, 0, 2*r)
def fill_ratio(r, h):
"""Horizontal pipe of radius r, filled with water at a height h.
Return the water fill ratio"""
c = chord(r, h-r)
theta = 2 * asin(c / (2*r))
factor = (theta - sin(theta)) / (2*pi)
return factor if h < r else 1-factor
# geometry
pd = 0.368 # m, pipe diameter
pl = 1.34 # m, pipe length
ao = 0.16 * 0.06 + pi/4 * 0.06**2 # m2, orifice area
# discharge coefficient
Cd = 0.8
# time traces
dt = 1e-2 # s
t = 0.0 # s
h = 1 * pd
time, height, fill = [t], [h], [fill_ratio(pd/2, h)]
while h > 1e-3:
if _debug: print(t)
t += dt
area = Cd * dt * (2*g*h)**0.5 * ao / pl
h -= calc_dh(area, h, pd/2)
time.append(t)
height.append(h)
fill.append(fill_ratio(pd/2, h))
plt.plot(time, fill)
plt.title('Spool Draining Time - 406 mm pipe')
plt.xlabel('Time [s]')
plt.ylabel('Fill ratio')
plt.xticks(range(11))
plt.yticks((0, 0.25, 0.5, 0.75, 1))
plt.grid()
plt.show()
print('%.1f seconds to empty the pipe' % t)