Georg_Brantegger/Python-DT_Slot_3
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Clean up before combining pipeline and

Browse Source 
reservoir code
This commit is contained in:
Georg ´Brantegger
2022-07-01 08:58:29 +02:00
parent d7585c3664
commit f774b2e52d
18 changed files with 0 additions and 8447 deletions
Download Patch File Download Diff File Expand all files Collapse all files

37
Ausgleichsbecken/dynamic_pipeline_pressure/Ausgleichsbecken.py
View File

@@ -1,31 +1,5 @@
import numpy as np
def Volume_trend(influx, outflux, timestep=1, V_0=0):
'''
Returns the trend and the volume and the final volume, defined
by influx and outflux patterns. The optional parameter timestep
defines the time increment over which the fluxes are changing.
'''
net_flux = influx-outflux
delta_V = net_flux*timestep
V_trend = V_0+np.cumsum(delta_V)
V_end = V_trend[-1]
return V_end, V_trend
def Height_trend(V_trend, area=1, h_crit_low=-np.inf, h_crit_high=np.inf):
'''
Returns the trend and the height and the final height, defined
by influx and outflux patterns as well as the crosssection area.
The optional parameters h_crit_low/high indicate limits that the height
should never exceed. If this occures, TRUE is returned in the corresponding
h_crit_flag.
'''
h_trend = V_trend/area
h_crit_flag_low = np.any(h_trend <= h_crit_low)
h_crit_flag_high = np.any(h_trend >= h_crit_high)
h_end = h_trend[-1]
return h_trend, h_end, h_crit_flag_low, h_crit_flag_high
def get_h_halfstep(initial_height, influx, outflux, timestep, area):
h0 = initial_height
Q_in = influx
@@ -45,7 +19,6 @@ def FODE_function(x, h, alpha, p, rho=1000., g=9.81):
def e_RK_4(yn, h, dt, Q0, Q1, A0, A1, p0, p1):
alpha = (A1/A0-1)
h_hs = get_h_halfstep(h, Q0, Q1, dt, A0)
p_hs = get_p_halfstep(p0, p1)
Y1 = yn
@@ -55,13 +28,3 @@ def e_RK_4(yn, h, dt, Q0, Q1, A0, A1, p0, p1):
ynp1 = yn + dt/6*(FODE_function(Y1, h, alpha, p)+2*FODE_function(Y2, h_hs, alpha, p_hs)+ \
2*FODE_function(Y3, h_hs, alpha, p_hs)+ FODE_function(Y4, h, alpha, p))
## testing
# if __name__ == "__main__":
# influx = np.full([1, 100], 6)
# outflux = np.full_like(influx, 4)
# V_end, V_trend = Volume_trend(influx, outflux, timestep=0.5, V_0 = 100)
# print(V_end)
# print(V_trend)

67
Ausgleichsbecken/static_pipeline_pressure/Ausgleichsbecken.py
View File

@@ -1,67 +0,0 @@
import numpy as np
def Volume_trend(influx, outflux, timestep=1, V_0=0):
'''
Returns the trend and the volume and the final volume, defined
by influx and outflux patterns. The optional parameter timestep
defines the time increment over which the fluxes are changing.
'''
net_flux = influx-outflux
delta_V = net_flux*timestep
V_trend = V_0+np.cumsum(delta_V)
V_end = V_trend[-1]
return V_end, V_trend
def Height_trend(V_trend, area=1, h_crit_low=-np.inf, h_crit_high=np.inf):
'''
Returns the trend and the height and the final height, defined
by influx and outflux patterns as well as the crosssection area.
The optional parameters h_crit_low/high indicate limits that the height
should never exceed. If this occures, TRUE is returned in the corresponding
h_crit_flag.
'''
h_trend = V_trend/area
h_crit_flag_low = np.any(h_trend <= h_crit_low)
h_crit_flag_high = np.any(h_trend >= h_crit_high)
h_end = h_trend[-1]
return h_trend, h_end, h_crit_flag_low, h_crit_flag_high
def get_h_halfstep(initial_height, influx, outflux, timestep, area):
h0 = initial_height
Q_in = influx
Q_out = outflux
dt = timestep
A = area
h_halfstep = h0+1/A*(Q_in-Q_out)*dt/2
def get_p_halfstep(p0, p1):
p_halfstep = (p0+p1)/2
def FODE_function(x, h, alpha, p, rho=1000., g=9.81):
f = x*abs(x)/h*alpha+g-p/(rho*h)
return f
def e_RK_4(yn, h, dt, Q0, Q1, A0, A1, p0, p1):
alpha = (A1/A0-1)
h_hs = get_h_halfstep(h, Q0, Q1, dt, A0)
p_hs = get_p_halfstep(p0, p1)
Y1 = yn
Y2 = yn + dt/2*FODE_function(Y1, h, alpha, p0)
Y3 = yn + dt/2*FODE_function(Y2, h_hs, alpha, p_hs)
Y4 = yn + dt*FODE_function(Y3, h_hs, alpha, p_hs)
ynp1 = yn + dt/6*(FODE_function(Y1, h, alpha, p)+2*FODE_function(Y2, h_hs, alpha, p_hs)+ \
2*FODE_function(Y3, h_hs, alpha, p_hs)+ FODE_function(Y4, h, alpha, p))
## testing
# if __name__ == "__main__":
# influx = np.full([1, 100], 6)
# outflux = np.full_like(influx, 4)
# V_end, V_trend = Volume_trend(influx, outflux, timestep=0.5, V_0 = 100)
# print(V_end)
# print(V_trend)

90
Ausgleichsbecken/static_pipeline_pressure/Ausgleichsbecken_class_file.py
View File

@@ -1,90 +0,0 @@
from Ausgleichsbecken import FODE_function, get_h_halfstep, get_p_halfstep
from pressure_conversion import pressure_conversion
class Ausgleichsbecken_class:
# units
area_unit = r'$\mathrm{m}^2$'
area_outflux_unit = r'$\mathrm{m}^2$'
level_unit = 'm'
volume_unit = r'$\mathrm{m}^3$'
flux_unit = r'$\mathrm{m}^3/\mathrm{s}$'
time_unit = 's'
pressure_unit = 'Pa'
# init
def __init__(self,area,outflux_area,level_min,level_max,timestep = 1):
self.area = area # base area of the rectangular structure
self.area_outflux = outflux_area # area of the outlet towards the pipeline
self.level_min = level_min # lowest allowed water level
self.level_max = level_max # highest allowed water level
self.timestep = timestep # timestep of the simulation
def update_volume(self):
self.volume = self.level*self.area
# setter
def set_initial_level(self,initial_level):
self.level = initial_level
self.update_volume()
def set_influx(self,influx):
self.influx = influx
def set_outflux(self,outflux):
self.outflux = outflux
# getter
def get_area(self):
print('The base area of the cuboid reservoir is', self.area, self.area_unit)
def get_outflux_area(self):
print('The outflux area from the cuboid reservoir to the pipeline is', \
self.area_outflux, self.area_outflux_unit)
def get_level(self):
print('The current level in the reservoir is', self.level , self.level_unit)
def get_crit_levels(self):
print('The critical water levels in the reservoir are: \n',\
' Minimum:', self.level_min , self.level_unit , '\n',\
' Maximum:', self.level_max , self.level_unit )
def get_volume(self):
print('The current water volume in the reservoir is', self.volume, self.volume_unit)
def get_timestep(self):
print('The timestep for the simulation is' , self.timestep, self.time_unit)
def get_influx(self):
print('The current influx is', self.influx, self.flux_unit)
def get_outflux(self):
print('The current outflux is', self.outflux, self.flux_unit)
# methods
def update_level(self,timestep):
# dont update volume here, because update_level gets called to calculate h_halfstep
net_flux = self.influx-self.outflux
delta_V = net_flux*timestep
new_level = (self.volume+delta_V)/self.area
return new_level
def e_RK_4(self):
# Update to deal with non constant pipeline pressure!
yn = self.outflux/self.area_outflux
h = self.level
dt = self.timestep
p,_ = pressure_conversion(self.initial_pressure,self.pressure_unit,'Pa')
p_hs,_ = pressure_conversion(self.initial_pressure,self.pressure_unit,'Pa')
alpha = (self.area_outflux/self.area-1)
h_hs = self.update_level(dt/2)
Y1 = yn
Y2 = yn + dt/2*FODE_function(Y1, h, alpha, self.initial_pressure)
Y3 = yn + dt/2*FODE_function(Y2, h_hs, alpha, p_hs)
Y4 = yn + dt*FODE_function(Y3, h_hs, alpha, p_hs)
ynp1 = yn + dt/6*(FODE_function(Y1, h, alpha, p)+2*FODE_function(Y2, h_hs, alpha, p_hs)+ \
2*FODE_function(Y3, h_hs, alpha, p_hs)+ FODE_function(Y4, h, alpha, p))
self.outflux = ynp1*self.area_outflux

145
Ausgleichsbecken/static_pipeline_pressure/Main_Program.ipynb
View File

@@ -1,145 +0,0 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from Ausgleichsbecken_class_file import Ausgleichsbecken_class\n",
"import matplotlib.pyplot as plt\n",
"from pressure_conversion import pressure_conversion"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# define constants\n",
"initial_level = 5. # m\n",
"initial_influx = 1. # m³/s\n",
"initial_outflux = 0. # m³/s\n",
"initial_pipeline_pressure = 1.\n",
"initial_pressure_unit = 'bar'\n",
"conversion_pressure_unit = 'mWS'\n",
"\n",
"area_base = 1. # m²\n",
"area_outflux = 0.5 # m²\n",
"critical_level_low = 0. # m\n",
"critical_level_high = 10. # m\n",
"simulation_timestep = 0.001 # s\n",
"\n",
"# for while loop\n",
"total_min_level = 0.01 # m\n",
"total_max_time = 150 # s"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib qt\n",
"\n",
"V = Ausgleichsbecken_class(area_base, area_outflux, critical_level_low, critical_level_high,simulation_timestep)\n",
"V.set_initial_level(initial_level) \n",
"V.set_influx(initial_influx)\n",
"V.set_outflux(initial_outflux)\n",
"\n",
"V.initial_pressure, V.pressure_unit = pressure_conversion(initial_pipeline_pressure,input_unit = initial_pressure_unit, target_unit = conversion_pressure_unit)\n",
"\n",
"time_vec = np.arange(0,total_max_time,simulation_timestep)\n",
"outflux_vec = np.empty_like(time_vec)\n",
"outflux_vec[0] = initial_outflux\n",
"level_vec = np.empty_like(time_vec)\n",
"level_vec[0] = initial_level\n",
" \n",
"i_max = -1\n",
"\n",
"for i in range(np.size(time_vec)-1):\n",
" V.e_RK_4()\n",
" V.level = V.update_level(V.timestep)\n",
" V.update_volume()\n",
" outflux_vec[i+1] = V.outflux\n",
" level_vec[i+1] = V.level\n",
" if V.level < total_min_level:\n",
" i_max = i\n",
" break\n",
"\n",
"\n",
"fig1, (ax1, ax2, ax3) = plt.subplots(3, 1)\n",
"fig1.set_figheight(10)\n",
"fig1.suptitle('Ausgleichsbecken')\n",
"\n",
"ax1.plot(time_vec[:i_max],level_vec[:i_max], label='Water level')\n",
"ax1.set_ylabel(r'$h$ ['+V.level_unit+']')\n",
"ax1.set_xlabel(r'$t$ ['+V.time_unit+']')\n",
"ax1.legend()\n",
"\n",
"ax2.plot(time_vec[:i_max],outflux_vec[:i_max], label='Outflux')\n",
"ax2.set_ylabel(r'$Q_{out}$ ['+V.flux_unit+']')\n",
"ax2.set_xlabel(r'$t$ ['+V.time_unit+']')\n",
"ax2.legend()\n",
"\n",
"# plt.subplots_adjust(left=0.2, bottom=0.2)\n",
"ax3.set_axis_off()\n",
"cell_text = np.array([[initial_level, V.level_unit], \\\n",
" [initial_influx, V.flux_unit], \\\n",
" [initial_outflux, V.flux_unit], \\\n",
" [round(V.initial_pressure,2), V.pressure_unit], \\\n",
" [simulation_timestep, V.time_unit], \\\n",
" [area_base, V.area_unit], \\\n",
" [area_outflux, V.area_unit]])\n",
"\n",
"row_labels =['initial_level', \\\n",
" 'initial_influx', \\\n",
" 'initial_outflux', \\\n",
" 'initial_pipeline_pressure', \\\n",
" 'simulation_timestep', \\\n",
" 'area_base', \\\n",
" 'area_outflux']\n",
"\n",
"plt.table(cellText=cell_text, \\\n",
" cellLoc='center', \\\n",
" colWidths=[0.3,0.1,0.3], \\\n",
" rowLabels=row_labels, \\\n",
" loc = 1, \\\n",
" rowLoc='left', \\\n",
" fontsize = 15.)\n",
"\n",
"fig1.tight_layout() "
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.8.13 ('Georg_DT_Slot3')",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.13"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "84fb123bdc47ab647d3782661abcbe80fbb79236dd2f8adf4cef30e8755eb2cd"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}

220
Ausgleichsbecken/static_pipeline_pressure/e-RK4-Test.ipynb
View File

File diff suppressed because one or more lines are too long

77
Ausgleichsbecken/static_pipeline_pressure/pressure_conversion.py
View File

@@ -1,77 +0,0 @@
# convert to Pa
def bar_to_pa(p):
return p*1e5
def mWS_to_pa(p):
return p*9.80665*1e3
def torr_to_pa(p):
return p*133.322
def atm_to_pa(p):
return p*101.325*1e3
def psi_to_pa(p):
return p*6894.8
# convert from Pa
def pa_to_bar(p):
return p*1e-5
def pa_to_mWS(p):
return p*1/(9.80665*1e3)
def pa_to_torr(p):
return p/133.322
def pa_to_atm(p):
return p*1/(101.325*1e3)
# converstion function
def pa_to_psi(p):
return p/6894.8
def pressure_conversion(pressure, input_unit = 'bar', target_unit = 'Pa'):
p = pressure
if input_unit.lower() == 'bar':
p_pa = bar_to_pa(p)
elif input_unit.lower() == 'mws':
p_pa = mWS_to_pa(p)
elif input_unit.lower() == 'torr':
p_pa = torr_to_pa(p)
elif input_unit.lower() == 'atm':
p_pa = atm_to_pa(p)
elif input_unit.lower() == 'psi':
p_pa = psi_to_pa(p)
elif input_unit.lower() == 'pa':
p_pa = p
else:
raise Exception('Given input unit not recognised. \n Known units are: Pa, bar, mWs, Torr, atm, psi')
if target_unit.lower() == 'bar':
return pa_to_bar(p_pa), target_unit
elif target_unit.lower() == 'mws':
return pa_to_mWS(p_pa), target_unit
elif target_unit.lower() == 'torr':
return pa_to_torr(p_pa), target_unit
elif target_unit.lower() == 'atm':
return pa_to_atm(p_pa), target_unit
elif target_unit.lower() =='psi':
return pa_to_psi(p_pa), target_unit
elif target_unit.lower() == 'pa':
return p_pa, target_unit
else:
raise Exception('Given target unit not recognised. \n Known units are: Pa, bar, mWs, Torr, atm, psi')
# testing_pressure_conversion
if __name__ == '__main__':
p = 1
unit_dict = ['Pa','Bar','Torr','Atm','MWS','psi']
for input_unit in unit_dict:
for target_unit in unit_dict:
converted_p = pressure_conversion(p,input_unit,target_unit)
print(input_unit,target_unit)
print(converted_p)