-
Notifications
You must be signed in to change notification settings - Fork 0
/
Thesis_ Severe_Slugging_Model.py
155 lines (120 loc) · 2.72 KB
/
Thesis_ Severe_Slugging_Model.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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
# -*- coding: utf-8 -*-
"""Untitled0.ipynb
Automatically generated by Colaboratory.
Original file is located at
https://colab.research.google.com/drive/1uCjl0klOEH5sc_tDZSAAsgb3mJ0VkqjE
"""
from math import *
import matplotlib.pyplot as plt
global p_res, Um_res, H_res, Nstep, T_res, z_res, H_increase, H_decrease, Usg_r_res, Usl_r_res
#Input - Flow inlet
Usgin=3
Uslin=0.3
#Input - Step
Nstep=50000
DeltaT=0.01
#Input - Geometry
L_l=100
D=0.045
A=pi*D**2/4
L_r=16.5
sinangle=1
#Input - Properties
p_out=100000
Rog=1
Rol=1000
visl=0.001
RT=100000
Pnormal=100000
#Input - Slip Relatiom
S=1.2
uo=0.35*sqrt(9.81*D)
#Input - Additional
valve_constant=0
#Initial value
p_res=[0]*Nstep
H_res=[0]*Nstep
z_res=[0]*Nstep
Um_res=[0]*Nstep
T_res=[0]*Nstep
H_increase=[0]*Nstep
H_decrease=[0]*Nstep
Usg_r_res=[0]*Nstep
Usl_r_res=[0]*Nstep
Rom=Rol
alfa_r=0
p=p_out+Rom*9.81*L_r
z=0.0001
U_level=0
Um=Uslin
H=1-alfa_r
for i in range(0,Nstep):
#Slug Length
L_s = -z+(1-alfa_r)*L_r
#Friction
Re=Rol*abs(Um)*D/visl
LambdaTurb=0.34*Re**(-0.25)
LambdaLam=64/Re
Lambda=max(LambdaTurb, LambdaLam)
Fric=0.5*Lambda*Rol*Um*abs(Um)*L_s/D
#Oscillation Friction
Lambda_stagnant=10
Fric_stagnant=0.5*Lambda_stagnant*Rol*U_level*abs(U_level)*L_s/D
#Valve Friction
Fric_valve=valve_constant*Rol*Um*abs(Um)
Fric=Fric+Fric_stagnant+Fric_valve
#Gravity
Grav_l=-Rol*9.81*abs(z*sinangle)
Grav_r=Rom*9.81*L_r
Grav=Grav_l+Grav_r
#Momentum Balance
Um_old=Um
Um=Um_old+DeltaT*((p-p_out)-Fric-Grav)/(-z*Rol+L_r*Rom)
#Slip Relation
Ug_r=S*Um+uo
#3 State Equation
Usg_l=(Um-Uslin)*(Um>0)*(z>0)
Usl_l=(Uslin*(z>0)+Um*(z<0))*(Um>0)+Um*(Um<0)
Usg_r=(alfa_r*Ug_r*(z>0)+Um*(z<0))*(Um>0)+Um*(Um<0)
Usl_r=(Um-alfa_r*Ug_r)*(Um>0)*(z>0)
U=Usg_l-Usg_r-Usl_l+Usl_r
#Alfa
alfa_r_old=alfa_r
alfa_r=alfa_r_old+0.5*DeltaT*U/(-z+L_r)
alfa_r=max(0,alfa_r)
alfa_r=min(1,alfa_r)
H_old=H
H=1-alfa_r
dHdt=(H-H_old)/DeltaT
if dHdt>0:
H_increase[i]=H
else:
H_decrease[i]=H
#Level
U_level=(-Uslin+Um)
U_level=U_level*(z<0)
z_old=z
z=z_old+DeltaT*U_level
z=(z<0)*z+(z>0)*0.001
#Pressure
p_old=p
p=p_old+DeltaT*(Usgin*Pnormal-Usg_l*p_old)/(L_l+z)-DeltaT*(p_old*U_level)/(L_l+z)
#Mixture Density
Rog=p/RT
Rom=alfa_r*Rog+(1-alfa_r)*Rol
#Results
p_res[i]=p
H_res[i]=(1-alfa_r)
Um_res[i]=Um
Usg_r_res[i]=Usg_r
Usl_r_res[i]=Usl_r
T_res[i]=i*DeltaT
z_res[i]=z
#Plot
x1 = T_res
y1 = Usg_r_res
fig, ax = plt.subplots()
ax.plot(x1, y1)
ax.set(xlabel='T', ylabel='Usl', title='Simplified Model')
ax.grid()
plt.show()