Merge branch 'Dev'
This commit is contained in:
Brantegger Georg
@@ -1,3 +1,4 @@
|
||||
from logging import exception
|
||||
import numpy as np
|
||||
|
||||
#importing pressure conversion function
|
||||
@@ -8,8 +9,19 @@ parent = os.path.dirname(current)
|
||||
sys.path.append(parent)
|
||||
from functions.pressure_conversion import pressure_conversion
|
||||
|
||||
def FODE_function(x, h, alpha, p, rho=1000., g=9.81):
|
||||
f = x*abs(x)/h*alpha+g-p/(rho*h)
|
||||
def FODE_function(x_out,h,A,A_a,p,rho,g):
|
||||
# (FODE ... first order differential equation)
|
||||
# based on the outflux formula by Andreas Malcherek
|
||||
# https://www.youtube.com/watch?v=8HO2LwqOhqQ
|
||||
# adapted for a pressurized pipeline into which the reservoir effuses
|
||||
# and flow direction
|
||||
# x_out ... effusion velocity
|
||||
# h ... level in the reservoir
|
||||
# A_a ... Outflux_Area
|
||||
# A ... Reservoir_Area
|
||||
# g ... gravitational acceleration
|
||||
# rho ... density of the liquid in the reservoir
|
||||
f = x_out*abs(x_out)/h*(A_a/A-1.)+g-p/(rho*h)
|
||||
return f
|
||||
|
||||
|
||||
@@ -19,71 +31,92 @@ class Ausgleichsbecken_class:
|
||||
# units are used to label graphs and print units are used to have a bearable format when using pythons print()
|
||||
area_unit = r'$\mathrm{m}^2$'
|
||||
area_outflux_unit = r'$\mathrm{m}^2$'
|
||||
density_unit = r'$\mathrm{kg}/\mathrm{m}^3$'
|
||||
flux_unit = r'$\mathrm{m}^3/\mathrm{s}$'
|
||||
level_unit = 'm'
|
||||
pressure_unit = 'Pa'
|
||||
time_unit = 's'
|
||||
velocity_unit = r'$\mathrm{m}/\mathrm{s}$'
|
||||
volume_unit = r'$\mathrm{m}^3$'
|
||||
|
||||
area_unit_print = 'm²'
|
||||
area_outflux_unit_print = 'm²'
|
||||
density_unit_print = 'kg/m³'
|
||||
flux_unit_print = 'm³/s'
|
||||
level_unit_print = 'm'
|
||||
pressure_unit_print = '--' # will be set by .set_pressure() method
|
||||
time_unit_print = 's'
|
||||
velocity_unit_print = 'm/s'
|
||||
volume_unit_print = 'm³'
|
||||
|
||||
g = 9.81 # m/s²
|
||||
rho = 1000 # kg/m³
|
||||
g = 9.81 # m/s² gravitational acceleration
|
||||
|
||||
|
||||
# init
|
||||
def __init__(self,area,outflux_area,level_min = 0,level_max = np.inf ,timestep = 1):
|
||||
def __init__(self,area,outflux_area,level_min = 0,level_max = np.inf ,timestep = 1,rho = 1000):
|
||||
self.area = area # base area of the rectangular structure
|
||||
self.area_outflux = outflux_area # area of the outlet towards the pipeline
|
||||
self.density = rho # density of the liquid in the system
|
||||
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
|
||||
|
||||
# initialize for get_info
|
||||
self.level = "--"
|
||||
self.influx = "--"
|
||||
self.level = "--"
|
||||
self.outflux = "--"
|
||||
self.volume = "--"
|
||||
|
||||
|
||||
# setter
|
||||
def set_volume(self):
|
||||
self.volume = self.level*self.area
|
||||
|
||||
def set_initial_level(self,initial_level):
|
||||
self.level = initial_level
|
||||
self.set_volume()
|
||||
# sets the level in the reservoir and should only be called during initialization
|
||||
if self.level == '--':
|
||||
self.level = initial_level
|
||||
else:
|
||||
raise Exception('Initial level was already set once. Use the .update_level(self,timestep) method to update level based on net flux.')
|
||||
|
||||
def set_level(self,level):
|
||||
self.level = level
|
||||
|
||||
def set_influx(self,influx):
|
||||
# sets influx to the reservoir in m³/s
|
||||
# positive influx means that liquid flows into the reservoir
|
||||
self.influx = influx
|
||||
|
||||
def set_outflux(self,outflux):
|
||||
self.outflux = outflux
|
||||
self.outflux_vel = outflux/self.area_outflux
|
||||
# sets outflux to the reservoir in m³/s
|
||||
# positive outflux means that liquid flows out of reservoir the reservoir
|
||||
self.outflux = outflux
|
||||
|
||||
def set_pressure(self,pressure,pressure_unit,display_pressure_unit):
|
||||
def set_initial_pressure(self,pressure,display_pressure_unit):
|
||||
# sets the static pressure present at the outlet of the reservoir
|
||||
# units are used to convert and display the pressure
|
||||
self.pressure = pressure
|
||||
self.pressure_unit = pressure_unit
|
||||
self.pressure_unit_print = display_pressure_unit
|
||||
|
||||
def set_steady_state(self,ss_influx,ss_level,pressure_unit,display_pressure_unit):
|
||||
ss_outflux = ss_influx
|
||||
ss_outflux_vel = ss_outflux/self.area_outflux
|
||||
ss_pressure = self.rho*self.g*ss_level-ss_outflux_vel**2*self.rho/2
|
||||
def set_pressure(self,pressure):
|
||||
# sets the static pressure present at the outlet of the reservoir
|
||||
self.pressure = pressure
|
||||
|
||||
def set_steady_state(self,ss_influx,ss_level,display_pressure_unit):
|
||||
# set the steady state (ss) condition in which the net flux is zero
|
||||
# set pressure acting on the outflux area so that the level stays constant
|
||||
ss_outflux = ss_influx
|
||||
ss_influx_vel = abs(ss_influx/self.area)
|
||||
ss_outflux_vel = abs(ss_outflux/self.area_outflux)
|
||||
ss_pressure = self.density*self.g*ss_level+self.density*ss_outflux_vel*(ss_influx_vel-ss_outflux_vel)
|
||||
|
||||
self.set_initial_level(ss_level)
|
||||
self.set_influx(ss_influx)
|
||||
self.set_initial_level(ss_level)
|
||||
self.set_initial_pressure(ss_pressure,display_pressure_unit)
|
||||
self.set_outflux(ss_outflux)
|
||||
self.set_pressure(ss_pressure,pressure_unit,display_pressure_unit)
|
||||
|
||||
# getter
|
||||
def get_info(self, full = False):
|
||||
new_line = '\n'
|
||||
p,_ = pressure_conversion(self.pressure,self.pressure_unit,self.pressure_unit_print)
|
||||
p = pressure_conversion(self.pressure,self.pressure_unit,self.pressure_unit_print)
|
||||
outflux_vel = self.outflux/self.area_outflux
|
||||
|
||||
|
||||
if full == True:
|
||||
@@ -98,9 +131,10 @@ class Ausgleichsbecken_class:
|
||||
f"Volume in reservoir = {self.volume:<10} {self.volume_unit_print} {new_line}"
|
||||
f"Current influx = {self.influx:<10} {self.flux_unit_print} {new_line}"
|
||||
f"Current outflux = {self.outflux:<10} {self.flux_unit_print} {new_line}"
|
||||
f"Current outflux vel = {round(self.outflux_vel,3):<10} {self.velocity_unit_print} {new_line}"
|
||||
f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_print} {new_line}"
|
||||
f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}"
|
||||
f"Simulation timestep = {self.timestep:<10} {self.time_unit_print} {new_line}"
|
||||
f"Density of liquid = {self.density:<10} {self.density_unit_print} {new_line}"
|
||||
f"----------------------------- {new_line}")
|
||||
else:
|
||||
# :<10 pads the self.value to be 10 characters wide
|
||||
@@ -110,38 +144,69 @@ class Ausgleichsbecken_class:
|
||||
f"Volume in reservoir = {self.volume:<10} {self.volume_unit_print} {new_line}"
|
||||
f"Current influx = {self.influx:<10} {self.flux_unit_print} {new_line}"
|
||||
f"Current outflux = {self.outflux:<10} {self.flux_unit_print} {new_line}"
|
||||
f"Current outflux vel = {round(self.outflux_vel,3):<10} {self.velocity_unit_print} {new_line}"
|
||||
f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_print} {new_line}"
|
||||
f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}"
|
||||
f"----------------------------- {new_line}")
|
||||
|
||||
print(print_str)
|
||||
|
||||
def get_current_level(self):
|
||||
return self.level
|
||||
|
||||
def get_current_influx(self):
|
||||
return self.influx
|
||||
|
||||
def get_current_outflux(self):
|
||||
return self.outflux
|
||||
|
||||
def get_current_pressure(self):
|
||||
return self.pressure
|
||||
|
||||
|
||||
|
||||
# methods
|
||||
|
||||
def update_level(self,timestep):
|
||||
# update level based on net flux and timestep by calculating the volume change in
|
||||
# the timestep and the converting the new volume to a level by assuming a cuboid reservoir
|
||||
|
||||
# cannot set new level directly in this method, because it gets called to calcuate during the Runge Kutta
|
||||
# to calculate a ficticious level at half the timestep
|
||||
net_flux = self.influx-self.outflux
|
||||
delta_V = net_flux*timestep
|
||||
new_level = (self.volume+delta_V)/self.area
|
||||
delta_level = net_flux*timestep/self.area
|
||||
new_level = (self.level+delta_level)
|
||||
return new_level
|
||||
|
||||
def update_pressure(self):
|
||||
influx_vel = abs(self.influx/self.area)
|
||||
outflux_vel = abs(self.outflux/self.area_outflux)
|
||||
p_new = self.density*self.g*self.level+self.density*outflux_vel*(influx_vel-outflux_vel)
|
||||
return p_new
|
||||
|
||||
def e_RK_4(self):
|
||||
def timestep_reservoir_evolution(self):
|
||||
# update outflux and outflux velocity based on current pipeline pressure and waterlevel in reservoir
|
||||
yn = self.outflux_vel
|
||||
h = self.level
|
||||
dt = self.timestep
|
||||
p = self.pressure
|
||||
# assume constant pipeline pressure during timestep (see comments in main_programm)
|
||||
p_hs = self.pressure
|
||||
alpha = (self.area_outflux/self.area-1)
|
||||
rho = self.density
|
||||
g = self.g
|
||||
A = self.area
|
||||
A_a = self.area_outflux
|
||||
yn = self.outflux/A_a # outflux velocity
|
||||
h = self.level
|
||||
h_hs = self.update_level(dt/2)
|
||||
p = self.pressure
|
||||
p_hs = self.pressure + rho*g*(h_hs-h)
|
||||
# explicit 4 step Runge Kutta
|
||||
Y1 = yn
|
||||
Y2 = yn + dt/2*FODE_function(Y1, h, alpha, self.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))
|
||||
Y2 = yn + dt/2*FODE_function(Y1,h,A,A_a,p,rho,g)
|
||||
Y3 = yn + dt/2*FODE_function(Y2,h_hs,A,A_a,p_hs,rho,g)
|
||||
Y4 = yn + dt*FODE_function(Y3,h_hs,A,A_a,p_hs,rho,g)
|
||||
ynp1 = yn + dt/6*(FODE_function(Y1,h,A,A_a,p,rho,g)+2*FODE_function(Y2,h_hs,A,A_a,p_hs,rho,g)+ \
|
||||
2*FODE_function(Y3,h_hs,A,A_a,p_hs,rho,g)+ FODE_function(Y4,h,A,A_a,p,rho,g))
|
||||
|
||||
self.outflux_vel = ynp1
|
||||
self.outflux = ynp1*self.area_outflux
|
||||
new_outflux = ynp1*A_a
|
||||
new_level = self.update_level(dt)
|
||||
new_pressure = self.update_pressure()
|
||||
|
||||
self.set_outflux(new_outflux)
|
||||
self.set_level(new_level)
|
||||
self.set_pressure(new_pressure)
|
||||
@@ -1,169 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"from Ausgleichsbecken_class_file import Ausgleichsbecken_class\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"\n",
|
||||
"#importing pressure conversion function\n",
|
||||
"import sys\n",
|
||||
"import os\n",
|
||||
"current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n",
|
||||
"parent = os.path.dirname(current)\n",
|
||||
"sys.path.append(parent)\n",
|
||||
"from functions.pressure_conversion import pressure_conversion"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# define constants\n",
|
||||
"initial_level = 10. # m\n",
|
||||
"initial_influx = 5. # m³/s\n",
|
||||
"initial_outflux = 1. # m³/s\n",
|
||||
"initial_pipeline_pressure = 10.\n",
|
||||
"initial_pressure_unit = 'mWS'\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 = 1000 # s"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"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",
|
||||
"# converted_pressure,_ = pressure_conversion(initial_pipeline_pressure,input_unit = initial_pressure_unit, target_unit = 'Pa')\n",
|
||||
"# V.pressure = converted_pressure\n",
|
||||
"V.set_steady_state(initial_influx,initial_level,initial_pressure_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] = V.outflux\n",
|
||||
"level_vec = np.empty_like(time_vec)\n",
|
||||
"level_vec[0] = V.level\n",
|
||||
"\n",
|
||||
"# pressure_vec = np.full_like(time_vec,converted_pressure)*((np.sin(time_vec)+1)*np.exp(-time_vec/50))\n",
|
||||
"pressure_vec = np.full_like(time_vec,V.pressure)\n",
|
||||
" \n",
|
||||
"i_max = -1\n",
|
||||
"\n",
|
||||
"for i in range(np.size(time_vec)-1):\n",
|
||||
" # update to include p_halfstep\n",
|
||||
" V.pressure = pressure_vec[i]\n",
|
||||
" V.e_RK_4()\n",
|
||||
" V.level = V.update_level(V.timestep)\n",
|
||||
" V.set_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"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\n",
|
||||
"fig1, (ax1, ax2, ax3, ax4) = plt.subplots(4, 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",
|
||||
"ax3.plot(time_vec[:i_max],pressure_conversion(pressure_vec[:i_max],'Pa',conversion_pressure_unit)[0], label='Pipeline pressure at reservoir')\n",
|
||||
"ax3.set_ylabel(r'$p_{pipeline}$ ['+conversion_pressure_unit+']')\n",
|
||||
"ax3.set_xlabel(r'$t$ ['+V.time_unit+']')\n",
|
||||
"ax3.legend()\n",
|
||||
"\n",
|
||||
"# plt.subplots_adjust(left=0.2, bottom=0.2)\n",
|
||||
"ax4.set_axis_off()\n",
|
||||
"cell_text = np.array([[level_vec[0], V.level_unit], \\\n",
|
||||
" [initial_influx, V.flux_unit], \\\n",
|
||||
" [outflux_vec[0], V.flux_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",
|
||||
" '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 ('DT_Slot_3')",
|
||||
"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": "4a28055eb8a3160fa4c7e4fca69770c4e0a1add985300856aa3fcf4ce32a2c48"
|
||||
}
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 2
|
||||
}
|
||||
198
Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb
Normal file
198
Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb
Normal file
@@ -0,0 +1,198 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"from Ausgleichsbecken_class_file import Ausgleichsbecken_class\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"\n",
|
||||
"#importing pressure conversion function\n",
|
||||
"import sys\n",
|
||||
"import os\n",
|
||||
"current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n",
|
||||
"parent = os.path.dirname(current)\n",
|
||||
"sys.path.append(parent)\n",
|
||||
"from functions.pressure_conversion import pressure_conversion"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"L = 1000.\n",
|
||||
"n = 10000 # number of pipe segments in discretization\n",
|
||||
"c = 400. \n",
|
||||
"dx = L/n # length of each pipe segment\n",
|
||||
"dt = dx/c \n",
|
||||
"\n",
|
||||
"# # define constants\n",
|
||||
"# initial_level = 10.1 # m\n",
|
||||
"# initial_influx = 0.8 # m³/s\n",
|
||||
"# conversion_pressure_unit = 'mWS'\n",
|
||||
"\n",
|
||||
"# area_base = 75. # m²\n",
|
||||
"# area_outflux = (0.9/2)**2*np.pi # m²\n",
|
||||
"# critical_level_low = 0. # m\n",
|
||||
"# critical_level_high = 10. # m\n",
|
||||
"# simulation_timestep = dt # s\n",
|
||||
"\n",
|
||||
"# # for while loop\n",
|
||||
"# total_min_level = 0.01 # m\n",
|
||||
"# total_max_time = 100 # s\n",
|
||||
"\n",
|
||||
"# nt = int(total_max_time//simulation_timestep)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# define constants\n",
|
||||
"initial_level = 10.1 # m\n",
|
||||
"initial_influx = 1. # m³/s\n",
|
||||
"# initial_outflux = 1. # m³/s\n",
|
||||
"# initial_pipeline_pressure = 10.\n",
|
||||
"# initial_pressure_unit = 'mWS'\n",
|
||||
"conversion_pressure_unit = 'mWS'\n",
|
||||
"\n",
|
||||
"area_base = 75. # m²\n",
|
||||
"area_outflux = 2. # m²\n",
|
||||
"critical_level_low = 0. # m\n",
|
||||
"critical_level_high = 10. # m\n",
|
||||
"simulation_timestep = dt # s\n",
|
||||
"\n",
|
||||
"# for while loop\n",
|
||||
"total_min_level = 0.01 # m\n",
|
||||
"total_max_time = 100 # s\n",
|
||||
"\n",
|
||||
"nt = int(total_max_time//simulation_timestep)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"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",
|
||||
"# V.set_initial_pressure(pressure_conversion(initial_pipeline_pressure,input_unit = initial_pressure_unit, target_unit = 'Pa'),conversion_pressure_unit)\n",
|
||||
"# V.pressure = converted_pressure\n",
|
||||
"V.set_steady_state(initial_influx,initial_level,conversion_pressure_unit)\n",
|
||||
"\n",
|
||||
"time_vec = np.arange(0,nt+1,1)*simulation_timestep\n",
|
||||
"outflux_vec = np.zeros_like(time_vec)\n",
|
||||
"outflux_vec[0] = V.get_current_outflux()\n",
|
||||
"level_vec = np.zeros_like(time_vec)\n",
|
||||
"level_vec[0] = V.get_current_level()\n",
|
||||
"pressure_vec = np.zeros_like(time_vec)\n",
|
||||
"pressure_vec[0] = V.get_current_pressure()\n",
|
||||
"\n",
|
||||
"# pressure_vec = np.full_like(time_vec,converted_pressure)*((np.sin(time_vec)+1)*np.exp(-time_vec/50))\n",
|
||||
" \n",
|
||||
"i_max = -1\n",
|
||||
"\n",
|
||||
"for i in range(1,nt+1):\n",
|
||||
" V.set_pressure(pressure_vec[i-1])\n",
|
||||
" V.set_outflux(outflux_vec[i-1])\n",
|
||||
" V.timestep_reservoir_evolution()\n",
|
||||
" outflux_vec[i] = V.get_current_outflux()\n",
|
||||
" level_vec[i] = V.get_current_level()\n",
|
||||
" pressure_vec[i] = V.get_current_pressure()\n",
|
||||
" if V.level < total_min_level:\n",
|
||||
" i_max = i\n",
|
||||
" break\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\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",
|
||||
"ax3.plot(time_vec[:i_max],pressure_conversion(pressure_vec[:i_max],'Pa',conversion_pressure_unit), label='Pipeline pressure at reservoir')\n",
|
||||
"ax3.set_ylabel(r'$p_{pipeline}$ ['+conversion_pressure_unit+']')\n",
|
||||
"ax3.set_xlabel(r'$t$ ['+V.time_unit+']')\n",
|
||||
"ax3.legend()\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"fig1.tight_layout() "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"10.1"
|
||||
]
|
||||
},
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"V.get_current_level()"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3.8.13 ('DT_Slot_3')",
|
||||
"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": "4a28055eb8a3160fa4c7e4fca69770c4e0a1add985300856aa3fcf4ce32a2c48"
|
||||
}
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 2
|
||||
}
|
||||
@@ -1,51 +1,42 @@
|
||||
import numpy as np
|
||||
|
||||
#importing pressure conversion function
|
||||
import sys
|
||||
import os
|
||||
current = os.path.dirname(os.path.realpath(__file__))
|
||||
parent = os.path.dirname(current)
|
||||
sys.path.append(parent)
|
||||
from functions.pressure_conversion import pressure_conversion
|
||||
|
||||
|
||||
|
||||
|
||||
class Druckrohrleitung_class:
|
||||
# units
|
||||
acceleration_unit = r'$\mathrm{m}/\mathrm{s}^2$'
|
||||
angle_unit = 'rad'
|
||||
area_unit = r'$\mathrm{m}^2$'
|
||||
density_unit = r'$\mathrm{kg}/\mathrm{m}^3$'
|
||||
flux_unit = r'$\mathrm{m}^3/\mathrm{s}$'
|
||||
length_unit = 'm'
|
||||
time_unit = 's'
|
||||
velocity_unit = r'$\mathrm{m}/\mathrm{s}$' # for flux and pressure propagation
|
||||
volume_unit = r'$\mathrm{m}^3$'
|
||||
acceleration_unit = r'$\mathrm{m}/\mathrm{s}^2$'
|
||||
angle_unit = 'rad'
|
||||
area_unit = r'$\mathrm{m}^2$'
|
||||
density_unit = r'$\mathrm{kg}/\mathrm{m}^3$'
|
||||
flux_unit = r'$\mathrm{m}^3/\mathrm{s}$'
|
||||
length_unit = 'm'
|
||||
pressure_unit = 'Pa'
|
||||
time_unit = 's'
|
||||
velocity_unit = r'$\mathrm{m}/\mathrm{s}$' # for flux and pressure propagation
|
||||
volume_unit = r'$\mathrm{m}^3$'
|
||||
|
||||
acceleration_unit_print = 'm/s²'
|
||||
angle_unit_print = 'rad'
|
||||
area_unit_print = 'm²'
|
||||
density_unit_print = 'kg/m³'
|
||||
flux_unit_print = 'm³/s'
|
||||
length_unit_print = 'm'
|
||||
time_unit_print = 's'
|
||||
velocity_unit_print = 'm/s' # for flux and pressure propagation
|
||||
volume_unit_print = 'm³'
|
||||
|
||||
acceleration_unit_print = 'm/s²'
|
||||
angle_unit_print = 'rad'
|
||||
area_unit_print = 'm²'
|
||||
density_unit_print = 'kg/m³'
|
||||
flux_unit_print = 'm³/s'
|
||||
length_unit_print = 'm'
|
||||
time_unit_print = 's'
|
||||
velocity_unit_print = 'm/s' # for flux and pressure propagation
|
||||
volume_unit_print = 'm³'
|
||||
|
||||
# init
|
||||
def __init__(self,total_length,diameter,number_segments,pipeline_angle,Darcy_friction_factor,rho=1000,g=9.81):
|
||||
self.length = total_length
|
||||
self.dia = diameter
|
||||
self.n_seg = number_segments
|
||||
self.angle = pipeline_angle
|
||||
self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient
|
||||
self.rho = rho
|
||||
self.g = g
|
||||
self.length = total_length # total length of the pipeline
|
||||
self.dia = diameter # diameter of the pipeline
|
||||
self.n_seg = number_segments # number of segments for the method of characteristics
|
||||
self.angle = pipeline_angle # angle of the pipeline
|
||||
self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient
|
||||
self.density = rho # density of the liquid in the pipeline
|
||||
self.g = g # gravitational acceleration
|
||||
|
||||
self.dx = total_length/number_segments
|
||||
self.l_vec = np.arange(0,(number_segments+1)*self.dx,self.dx)
|
||||
self.A = (diameter/2)**2*np.pi
|
||||
|
||||
self.dx = total_length/number_segments # length of each segment
|
||||
self.l_vec = np.arange(0,(number_segments+1),1)*self.dx # vector giving the distance from each node to the start of the pipeline
|
||||
|
||||
# initialize for get_info method
|
||||
self.c = '--'
|
||||
@@ -53,32 +44,32 @@ class Druckrohrleitung_class:
|
||||
|
||||
# setter
|
||||
def set_pressure_propagation_velocity(self,c):
|
||||
self.c = c
|
||||
self.dt = self.dx/c
|
||||
self.c = c # propagation velocity of the pressure wave
|
||||
self.dt = self.dx/c # timestep derived from c, demanded by the method of characteristics
|
||||
|
||||
def set_number_of_timesteps(self,number_timesteps):
|
||||
self.nt = number_timesteps
|
||||
self.nt = number_timesteps # number of timesteps
|
||||
if self.c == '--':
|
||||
raise Exception('Please set the pressure propagation velocity before setting the number of timesteps.')
|
||||
else:
|
||||
self.t_vec = np.arange(0,self.nt*self.dt,self.dt)
|
||||
|
||||
def set_initial_pressure(self,pressure,pressure_unit,display_pressure_unit):
|
||||
if np.size(pressure) == 1:
|
||||
def set_initial_pressure(self,pressure):
|
||||
# initialize the pressure distribution in the pipeline
|
||||
if np.size(pressure) == 1:
|
||||
self.p0 = np.full_like(self.l_vec,pressure)
|
||||
elif np.size(pressure) == np.size(self.l_vec):
|
||||
self.p0 = pressure
|
||||
else:
|
||||
raise Exception('Unable to assign initial pressure. Input has to be of size 1 or' + np.size(self.l_vec))
|
||||
self.pressure_unit = pressure_unit
|
||||
self.pressure_unit_print = display_pressure_unit
|
||||
|
||||
#initialize the vectors in which the old and new pressures are stored for the method of characteristics
|
||||
self.p_old = self.p0.copy()
|
||||
self.p = np.empty_like(self.p_old)
|
||||
self.p_old = self.p0.copy()
|
||||
self.p = self.p0.copy()
|
||||
|
||||
def set_initial_flow_velocity(self,velocity):
|
||||
if np.size(velocity) == 1:
|
||||
# initialize the velocity distribution in the pipeline
|
||||
if np.size(velocity) == 1:
|
||||
self.v0 = np.full_like(self.l_vec,velocity)
|
||||
elif np.size(velocity) == np.size(self.l_vec):
|
||||
self.v0 = velocity
|
||||
@@ -86,35 +77,54 @@ class Druckrohrleitung_class:
|
||||
raise Exception('Unable to assign initial velocity. Input has to be of size 1 or' + np.size(self.l_vec))
|
||||
|
||||
#initialize the vectors in which the old and new velocities are stored for the method of characteristics
|
||||
self.v_old = self.v0.copy()
|
||||
self.v = np.empty_like(self.v_old)
|
||||
self.v_old = self.v0.copy()
|
||||
self.v = self.v0.copy()
|
||||
|
||||
def set_boundary_conditions_next_timestep(self,v_reservoir,p_reservoir,v_turbine):
|
||||
rho = self.rho
|
||||
c = self.c
|
||||
f_D = self.f_D
|
||||
dt = self.dt
|
||||
D = self.dia
|
||||
g = self.g
|
||||
alpha = self.angle
|
||||
p_old = self.p_old[-2] # @ second to last node (the one before the turbine)
|
||||
v_old = self.v_old[-2] # @ second to last node (the one before the turbine)
|
||||
self.v_boundary_res = v_reservoir # at new timestep
|
||||
self.v_boundary_tur = v_turbine # at new timestep
|
||||
self.p_boundary_res = p_reservoir
|
||||
self.p_boundary_tur = p_old-rho*c*(v_turbine-v_old)+rho*c*dt*g*np.sin(alpha)-f_D*rho*c*dt/(2*D)*abs(v_old)*v_old
|
||||
self.v[0] = self.v_boundary_res.copy()
|
||||
self.v[-1] = self.v_boundary_tur.copy()
|
||||
self.p[0] = self.p_boundary_res.copy()
|
||||
self.p[-1] = self.p_boundary_tur.copy()
|
||||
def set_boundary_conditions_next_timestep(self,p_reservoir,v_turbine):
|
||||
# derived from the method of characteristics, one can set the boundary conditions for the pressures and flow velocities at the reservoir and the turbine
|
||||
# the boundary velocity at the turbine is specified by the flux through the turbine or an external boundary condition
|
||||
# the pressure at the turbine will be calculated using the forward characteristic
|
||||
# the boundary pressure at the reservoir is specified by the level in the reservoir of an external boundary condition
|
||||
# the velocity at the reservoir will be calculated using the backward characteristic
|
||||
|
||||
# constants for a cleaner formula
|
||||
rho = self.density
|
||||
c = self.c
|
||||
f_D = self.f_D
|
||||
dt = self.dt
|
||||
D = self.dia
|
||||
g = self.g
|
||||
alpha = self.angle
|
||||
p_old_tur = self.p_old[-2] # @ second to last node (the one before the turbine)
|
||||
v_old_tur = self.v_old[-2] # @ second to last node (the one before the turbine)
|
||||
p_old_res = self.p_old[1] # @ second node (the one after the reservoir)
|
||||
v_old_res = self.v_old[1] # @ second node (the one after the reservoir)
|
||||
# set the boundary conditions derived from reservoir and turbine
|
||||
v_boundary_tur = v_turbine # at new timestep
|
||||
p_boundary_res = p_reservoir # at new timestep
|
||||
# calculate the missing boundary conditions
|
||||
v_boundary_res = v_old_res+1/(rho*c)*(p_boundary_res-p_old_res)+dt*g*np.sin(alpha)-f_D*dt/(2*D)*abs(v_old_res)*v_old_res
|
||||
p_boundary_tur = p_old_tur-rho*c*(v_boundary_tur-v_old_tur)+rho*c*dt*g*np.sin(alpha)-f_D*rho*c*dt/(2*D)*abs(v_old_tur)*v_old_tur
|
||||
|
||||
def set_steady_state(self,ss_flux,ss_level_reservoir,pl_vec,h_vec,pressure_unit,display_pressure_unit):
|
||||
ss_v0 = np.full(self.n_seg+1,ss_flux/(self.dia**2/4*np.pi))
|
||||
ss_pressure = (self.rho*self.g*(ss_level_reservoir+h_vec)-ss_v0**2*self.rho/2)-(self.f_D*pl_vec/self.dia*self.rho/2*ss_v0**2)
|
||||
# write boundary conditions to the velocity/pressure vectors of the next timestep
|
||||
self.v[0] = v_boundary_res
|
||||
self.v[-1] = v_boundary_tur
|
||||
self.p[0] = p_boundary_res
|
||||
self.p[-1] = p_boundary_tur
|
||||
|
||||
def set_steady_state(self,ss_flux,ss_level_reservoir,area_reservoir,pl_vec,h_vec):
|
||||
# set the pressure and velocity distributions, that allow a constant flow of water from the (steady-state) reservoir to the (steady-state) turbine
|
||||
# the flow velocity is given by the constant flow through the pipe
|
||||
ss_v0 = np.full(self.n_seg+1,ss_flux/self.A)
|
||||
|
||||
# the static pressure is given by static state pressure of the reservoir, corrected for the hydraulic head of the pipe and friction losses
|
||||
ss_v_in_res = abs(ss_flux/area_reservoir)
|
||||
ss_v_out_res = abs(ss_flux/self.A)
|
||||
ss_pressure_res = self.density*self.g*(ss_level_reservoir)+self.density*ss_v_out_res*(ss_v_in_res-ss_v_out_res)
|
||||
ss_pressure = ss_pressure_res+(self.density*self.g*h_vec)-(self.f_D*pl_vec/self.dia*self.density/2*ss_v0**2)
|
||||
|
||||
self.set_initial_flow_velocity(ss_v0)
|
||||
self.set_initial_pressure(ss_pressure,pressure_unit,display_pressure_unit)
|
||||
self.set_initial_pressure(ss_pressure)
|
||||
|
||||
# getter
|
||||
def get_info(self):
|
||||
@@ -142,30 +152,27 @@ class Druckrohrleitung_class:
|
||||
f"Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object")
|
||||
|
||||
print(print_str)
|
||||
|
||||
|
||||
def get_boundary_conditions_next_timestep(self):
|
||||
print('The pressure at the reservoir for the next timestep is', '\n', \
|
||||
pressure_conversion(self.p_boundary_res,self.pressure_unit,self.pressure_unit_print), '\n', \
|
||||
'The velocity at the reservoir for the next timestep is', '\n', \
|
||||
self.v_boundary_res, self.velocity_unit_print, '\n', \
|
||||
'The pressure at the turbine for the next timestep is', '\n', \
|
||||
pressure_conversion(self.p_boundary_tur,self.pressure_unit,self.pressure_unit_print), '\n', \
|
||||
'The velocity at the turbine for the next timestep is', '\n', \
|
||||
self.v_boundary_tur, self.velocity_unit_print)
|
||||
def get_current_pressure_distribution(self):
|
||||
return self.p
|
||||
|
||||
def get_current_velocity_distribution(self):
|
||||
return self.v
|
||||
|
||||
def timestep_characteristic_method(self):
|
||||
#number of nodes
|
||||
nn = self.n_seg+1
|
||||
rho = self.rho
|
||||
c = self.c
|
||||
f_D = self.f_D
|
||||
dt = self.dt
|
||||
D = self.dia
|
||||
g = self.g
|
||||
alpha = self.angle
|
||||
# use the method of characteristics to calculate the pressure and velocities at all nodes except the boundary ones
|
||||
# they are set with the .set_boundary_conditions_next_timestep() method beforehand
|
||||
|
||||
nn = self.n_seg+1 # number of nodes
|
||||
rho = self.density # density of liquid
|
||||
c = self.c # pressure propagation velocity
|
||||
f_D = self.f_D # Darcy friction coefficient
|
||||
dt = self.dt # timestep
|
||||
D = self.dia # pipeline diametet
|
||||
g = self.g # graviational acceleration
|
||||
alpha = self.angle # pipeline angle
|
||||
|
||||
# Vectorize this loop?
|
||||
for i in range(1,nn-1):
|
||||
self.v[i] = 0.5*(self.v_old[i+1]+self.v_old[i-1])-0.5/(rho*c)*(self.p_old[i+1]-self.p_old[i-1]) \
|
||||
+dt*g*np.sin(alpha)-f_D*dt/(4*D)*(abs(self.v_old[i+1])*self.v_old[i+1]+abs(self.v_old[i-1])*self.v_old[i-1])
|
||||
@@ -173,5 +180,8 @@ class Druckrohrleitung_class:
|
||||
self.p[i] = 0.5*(self.p_old[i+1]+self.p_old[i-1]) - 0.5*rho*c*(self.v_old[i+1]-self.v_old[i-1]) \
|
||||
+f_D*rho*c*dt/(4*D)*(abs(self.v_old[i+1])*self.v_old[i+1]-abs(self.v_old[i-1])*self.v_old[i-1])
|
||||
|
||||
# prepare for next call
|
||||
# use .copy() to write data to another memory location and avoid the usual python reference pointer
|
||||
# else one can overwrite data by accidient and change two variables at once without noticing
|
||||
self.p_old = self.p.copy()
|
||||
self.v_old = self.v.copy()
|
||||
|
||||
@@ -1,237 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"from Druckrohrleitung_class_file import Druckrohrleitung_class\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"\n",
|
||||
"#importing pressure conversion function\n",
|
||||
"import sys\n",
|
||||
"import os\n",
|
||||
"current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n",
|
||||
"parent = os.path.dirname(current)\n",
|
||||
"sys.path.append(parent)\n",
|
||||
"from functions.pressure_conversion import pressure_conversion"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%matplotlib qt5\n",
|
||||
"#define constants\n",
|
||||
"\n",
|
||||
"g = 9.81 # gravitational acceleration [m/s²]\n",
|
||||
"rho = 1000 # density of water [kg/m³]\n",
|
||||
"\n",
|
||||
"L = 1000 # length of pipeline [m]\n",
|
||||
"D = 1 # pipe diameter [m]\n",
|
||||
"Q0 = 2 # initial flow in whole pipe [m³/s]\n",
|
||||
"h_res = 20 # water level in upstream reservoir [m]\n",
|
||||
"n = 10 # number of pipe segments in discretization\n",
|
||||
"nt = 100 # number of time steps after initial conditions\n",
|
||||
"f_D = 0.01 # Darcy friction factor\n",
|
||||
"c = 400 # propagation velocity of the pressure wave [m/s]\n",
|
||||
"h_pipe = 200 # hydraulic head without reservoir [m] \n",
|
||||
"alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# preparing the discretization and initial conditions\n",
|
||||
"\n",
|
||||
"dx = L/n # length of each pipe segment\n",
|
||||
"dt = dx/c # timestep according to method of characterisitics\n",
|
||||
"nn = n+1 # number of nodes\n",
|
||||
"pl_vec = np.arange(0,nn*dx,dx) # pl = pipe-length. position of the nodes on the pipeline\n",
|
||||
"t_vec = np.arange(0,nt*dt,dt) # time vector\n",
|
||||
"h_vec = np.arange(0,h_pipe+h_pipe/n,h_pipe/n) # hydraulic head of pipeline at each node\n",
|
||||
"\n",
|
||||
"v_init = np.full(nn,Q0/(D**2/4*np.pi))\n",
|
||||
"p_init = (rho*g*(h_res+h_vec)-v_init**2*rho/2)-(f_D*pl_vec/D*rho/2*v_init**2) # ref Wikipedia: Darcy Weisbach\n",
|
||||
"\n",
|
||||
"# storage vectors for old parameters\n",
|
||||
"v_old = v_init.copy()\n",
|
||||
"p_old = p_init.copy() \n",
|
||||
"\n",
|
||||
"# storage vectors for new parameters\n",
|
||||
"v_new = np.empty_like(v_old)\n",
|
||||
"p_new = np.empty_like(p_old)\n",
|
||||
"\n",
|
||||
"# storage vector for time evolution of parameters at node 0 (at reservoir)\n",
|
||||
"p_0 = np.full_like(t_vec,p_init[0])\n",
|
||||
"v_0 = np.full_like(t_vec,v_init[0])\n",
|
||||
"\n",
|
||||
"# storage vector for time evolution of parameters at node N+1 (at valve)\n",
|
||||
"p_np1 = np.full_like(t_vec,p_init[-1])\n",
|
||||
"v_np1 = np.full_like(t_vec,v_init[-1])\n",
|
||||
"\n",
|
||||
"for it in range(1,nt):\n",
|
||||
"\n",
|
||||
" # set boundary conditions\n",
|
||||
" v_new[-1] = 0 # in front of the instantaneously closing valve, the velocity is 0\n",
|
||||
" p_new[0] = p_init[0] # hydrostatic pressure from the reservoir\n",
|
||||
"\n",
|
||||
" # calculate the new parameters at first and last node\n",
|
||||
" v_new[0] = v_old[1]+1/(rho*c)*(p_init[0]-p_old[1])+dt*g*np.sin(alpha)-f_D*dt/(2*D)*abs(v_old[1])*v_old[1]\n",
|
||||
" p_new[-1] = p_old[-2]+rho*c*v_old[-2]-rho*c*f_D*dt/(2*D) *abs(v_old[-2])*v_old[-2]\n",
|
||||
"\n",
|
||||
" # calculate parameters at second to second-to-last nodes \n",
|
||||
" #equation 2-30 plus 2-31 (and refactor for v_i^j+1) in block 2\n",
|
||||
"\n",
|
||||
" for i in range(1,nn-1):\n",
|
||||
" v_new[i] = 0.5*(v_old[i-1]+v_old[i+1])+0.5/(rho*c)*(p_old[i-1]-p_old[i+1]) \\\n",
|
||||
" +dt*g*np.sin(alpha)-f_D*dt/(4*D)*(abs(v_old[i-1])*v_old[i-1]+abs(v_old[i+1])*v_old[i+1])\n",
|
||||
"\n",
|
||||
" p_new[i] = 0.5*rho*c*(v_old[i-1]-v_old[i+1])+0.5*(p_old[i-1]+p_old[i+1]) \\\n",
|
||||
" -rho*c*f_D*dt/(4*D)*(abs(v_old[i-1])*v_old[i-1]-abs(v_old[i+1])*v_old[i+1])\n",
|
||||
" \n",
|
||||
"\n",
|
||||
" # prepare for next loop\n",
|
||||
" # use .copy() to avoid that memory address is overwritten and hell breaks loose :D\n",
|
||||
" #https://www.geeksforgeeks.org/array-copying-in-python/\n",
|
||||
" p_old = p_new.copy()\n",
|
||||
" v_old = v_new.copy()\n",
|
||||
"\n",
|
||||
" # store parameters of node 1 (at reservoir)\n",
|
||||
" p_0[it] = p_new[0]\n",
|
||||
" v_0[it] = v_new[0]\n",
|
||||
" # store parameters of node N+1 (at reservoir)\n",
|
||||
" p_np1[it] = p_new[-1]\n",
|
||||
" v_np1[it] = v_new[-1]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n",
|
||||
"\n",
|
||||
"pipe.set_pressure_propagation_velocity(c)\n",
|
||||
"pipe.set_number_of_timesteps(nt)\n",
|
||||
"\n",
|
||||
"pipe.set_initial_pressure(p_init)\n",
|
||||
"pipe.set_initial_flow_velocity(v_init)\n",
|
||||
"pipe.set_boundary_conditions_next_timestep(v_0[0],p_0[0],v_np1[0])\n",
|
||||
"\n",
|
||||
"# storage vector for time evolution of parameters at node 0 (at reservoir)\n",
|
||||
"pipe.p_0 = np.full_like(t_vec,p_init[0])\n",
|
||||
"pipe.v_0 = np.full_like(t_vec,v_init[0])\n",
|
||||
"\n",
|
||||
"# storage vector for time evolution of parameters at node N+1 (at valve)\n",
|
||||
"pipe.p_np1 = np.full_like(t_vec,p_init[-1])\n",
|
||||
"pipe.v_np1 = np.full_like(t_vec,v_init[-1])\n",
|
||||
"\n",
|
||||
"fig2,axs2 = plt.subplots(2,1)\n",
|
||||
"axs2[0].set_title('Pressure distribution in pipeline')\n",
|
||||
"axs2[1].set_title('Velocity distribution in pipeline')\n",
|
||||
"axs2[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
|
||||
"axs2[0].set_ylabel(r'$p$ [mWS]')\n",
|
||||
"axs2[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
|
||||
"axs2[1].set_ylabel(r'$p$ [mWS]')\n",
|
||||
"lo_00, = axs2[0].plot(pl_vec,pressure_conversion(pipe.p_old,'Pa','mWS')[0],marker='.')\n",
|
||||
"lo_01, = axs2[1].plot(pl_vec,pipe.v_old,marker='.')\n",
|
||||
"axs2[0].set_ylim([-2*np.max(pressure_conversion(p_init,'Pa','mWS')[0]),2*np.max(pressure_conversion(p_init,'Pa','mWS')[0])])\n",
|
||||
"axs2[1].set_ylim([-2*np.max(v_init),2*np.max(v_init)])\n",
|
||||
"fig2.tight_layout()\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"for it in range(1,pipe.nt):\n",
|
||||
" pipe.set_boundary_conditions_next_timestep(v_0[it],p_0[it],v_np1[it])\n",
|
||||
" pipe.timestep_characteristic_method()\n",
|
||||
" lo_00.set_ydata(pressure_conversion(pipe.p,'Pa','mWS')[0])\n",
|
||||
" lo_01.set_ydata(pipe.v)\n",
|
||||
"\n",
|
||||
" # store parameters of node 0 (at reservoir)\n",
|
||||
" pipe.p_0[it] = pipe.p[0]\n",
|
||||
" pipe.v_0[it] = pipe.v[0]\n",
|
||||
" # store parameters of node N+1 (at reservoir)\n",
|
||||
" pipe.p_np1[it] = pipe.p[-1]\n",
|
||||
" pipe.v_np1[it] = pipe.v[-1]\n",
|
||||
" \n",
|
||||
" fig2.suptitle(str(it))\n",
|
||||
" fig2.canvas.draw()\n",
|
||||
" fig2.tight_layout()\n",
|
||||
" plt.pause(0.2)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"fig3,axs3 = plt.subplots(2,2)\n",
|
||||
"axs3[0,0].plot(t_vec,pressure_conversion(pipe.p_0,'Pa','mWS')[0])\n",
|
||||
"axs3[0,1].plot(t_vec,pipe.v_0)\n",
|
||||
"axs3[1,0].plot(t_vec,pressure_conversion(pipe.p_np1,'Pa','mWS')[0])\n",
|
||||
"axs3[1,1].plot(t_vec,pipe.v_np1)\n",
|
||||
"axs3[0,0].set_title('Pressure Reservoir')\n",
|
||||
"axs3[0,1].set_title('Velocity Reservoir')\n",
|
||||
"axs3[1,0].set_title('Pressure Turbine')\n",
|
||||
"axs3[1,1].set_title('Velocity Turbine')\n",
|
||||
"axs3[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs3[0,0].set_ylabel(r'$p$ [mWS]')\n",
|
||||
"axs3[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs3[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
|
||||
"axs3[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs3[1,0].set_ylabel(r'$p$ [mWS]')\n",
|
||||
"axs3[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs3[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
|
||||
"fig3.tight_layout()\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"0.29590621205048523\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"print(np.mean(v_0))"
|
||||
]
|
||||
}
|
||||
],
|
||||
"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
|
||||
}
|
||||
278
Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb
Normal file
278
Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb
Normal file
@@ -0,0 +1,278 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"from Druckrohrleitung_class_file import Druckrohrleitung_class\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"\n",
|
||||
"#importing pressure conversion function\n",
|
||||
"import sys\n",
|
||||
"import os\n",
|
||||
"current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n",
|
||||
"parent = os.path.dirname(current)\n",
|
||||
"sys.path.append(parent)\n",
|
||||
"from functions.pressure_conversion import pressure_conversion\n",
|
||||
"from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%matplotlib qt5\n",
|
||||
"#define constants pipe\n",
|
||||
"\n",
|
||||
"g = 9.81 # gravitational acceleration [m/s²]\n",
|
||||
"rho = 1000. # density of water [kg/m³]\n",
|
||||
"\n",
|
||||
"L = 1000. # length of pipeline [m]\n",
|
||||
"D = 0.9 # pipe diameter [m]\n",
|
||||
"h_res = 10. # water level in upstream reservoir [m]\n",
|
||||
"n = 50 # number of pipe segments in discretization\n",
|
||||
"nt = 1000 # number of time steps after initial conditions\n",
|
||||
"f_D = 0.01 # Darcy friction factor\n",
|
||||
"c = 400. # propagation velocity of the pressure wave [m/s]\n",
|
||||
"h_pipe = 105. # hydraulic head without reservoir [m] \n",
|
||||
"alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# preparing the discretization and initial conditions\n",
|
||||
"initial_flux = 0.8 # m³/s\n",
|
||||
"initial_level = h_res # m\n",
|
||||
"dx = L/n # length of each pipe segment\n",
|
||||
"dt = dx/c # timestep according to method of characterisitics\n",
|
||||
"nn = n+1 # number of nodes\n",
|
||||
"pl_vec = np.arange(0,nn,1)*dx # pl = pipe-length. position of the nodes on the pipeline\n",
|
||||
"t_vec = np.arange(0,nt,1)*dt # time vector\n",
|
||||
"h_vec = np.arange(0,nn,1)*h_pipe/n # hydraulic head of pipeline at each node\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# define constants reservoir\n",
|
||||
"conversion_pressure_unit = 'mWS'\n",
|
||||
"\n",
|
||||
"area_base = 75. # m²\n",
|
||||
"area_pipe = (D/2)**2*np.pi # m²\n",
|
||||
"critical_level_low = 0. # m\n",
|
||||
"critical_level_high = 100. # m\n",
|
||||
"\n",
|
||||
"# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n",
|
||||
"nt_eRK4 = 1 # number of simulation steps of reservoir in between timesteps of pipeline \n",
|
||||
"simulation_timestep = dt/nt_eRK4"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"V = Ausgleichsbecken_class(area_base,area_pipe,critical_level_low,critical_level_high,simulation_timestep)\n",
|
||||
"V.set_steady_state(initial_flux,initial_level,conversion_pressure_unit)\n",
|
||||
"\n",
|
||||
"pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n",
|
||||
"pipe.set_pressure_propagation_velocity(c)\n",
|
||||
"pipe.set_number_of_timesteps(nt)\n",
|
||||
"pipe.set_steady_state(initial_flux,initial_level,area_base,pl_vec,h_vec)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"The current attributes are: \n",
|
||||
"----------------------------- \n",
|
||||
"Current level = 10.0 m\n",
|
||||
"Volume in reservoir = -- m³ \n",
|
||||
"Current influx = 0.8 m³/s \n",
|
||||
"Current outflux = 0.8 m³/s \n",
|
||||
"Current outflux vel = 1.258 m/s \n",
|
||||
"Current pipe pressure = 9.844 mWS \n",
|
||||
"----------------------------- \n",
|
||||
"\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"V.get_info()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# initialization for timeloop\n",
|
||||
"\n",
|
||||
"level_vec = np.zeros_like(t_vec)\n",
|
||||
"level_vec[0] = V.get_current_level()\n",
|
||||
"\n",
|
||||
"# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n",
|
||||
"v_old = pipe.get_current_velocity_distribution()\n",
|
||||
"p_old = pipe.get_current_pressure_distribution()\n",
|
||||
"\n",
|
||||
"# prepare the vectors in which the temporal evolution of the boundary conditions are stored\n",
|
||||
" # keep in mind, that the velocity at the turbine and the pressure at the reservoir are set manually and\n",
|
||||
" # through the time evolution of the reservoir respectively \n",
|
||||
" # the pressure at the turbine and the velocity at the reservoir are calculated from the method of characteristics\n",
|
||||
"v_boundary_res = np.zeros_like(t_vec)\n",
|
||||
"v_boundary_tur = np.zeros_like(t_vec)\n",
|
||||
"p_boundary_res = np.zeros_like(t_vec)\n",
|
||||
"p_boundary_tur = np.zeros_like(t_vec)\n",
|
||||
"\n",
|
||||
"# set the boundary conditions for the first timestep\n",
|
||||
"v_boundary_res[0] = v_old[0]\n",
|
||||
"v_boundary_tur[0] = v_old[-1] \n",
|
||||
"p_boundary_res[0] = p_old[0]\n",
|
||||
"p_boundary_tur[0] = p_old[-1]\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"fig1,axs1 = plt.subplots(2,1)\n",
|
||||
"axs1[0].set_title('Pressure distribution in pipeline')\n",
|
||||
"axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
|
||||
"axs1[0].set_ylabel(r'$p$ [mWS]')\n",
|
||||
"lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,'Pa',conversion_pressure_unit),marker='.')\n",
|
||||
"axs1[0].set_ylim([0.9*np.min(pressure_conversion(p_old,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_old,'Pa',conversion_pressure_unit))])\n",
|
||||
"\n",
|
||||
"axs1[1].set_title('Velocity distribution in pipeline')\n",
|
||||
"axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
|
||||
"axs1[1].set_ylabel(r'$v$ [m/s]')\n",
|
||||
"lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n",
|
||||
"# axs1[1].set_ylim([0.9*np.min(v_old),1.1*np.max(v_boundary_res)])\n",
|
||||
"\n",
|
||||
"fig1.tight_layout()\n",
|
||||
"plt.pause(1)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\n",
|
||||
"for it_pipe in range(1,nt):\n",
|
||||
"# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n",
|
||||
" # set initial conditions for the reservoir time evolution calculted with e-RK4\n",
|
||||
" V.set_pressure(p_old[0])\n",
|
||||
" V.set_outflux(v_old[0]*area_pipe)\n",
|
||||
" # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n",
|
||||
" for it_res in range(nt_eRK4):\n",
|
||||
" V.timestep_reservoir_evolution() \n",
|
||||
" level_vec[it_pipe] = V.get_current_level() \n",
|
||||
"\n",
|
||||
" \n",
|
||||
" # set boundary conditions for the next timestep of the characteristic method\n",
|
||||
" p_boundary_res[it_pipe] = V.get_current_pressure()\n",
|
||||
" v_boundary_tur[it_pipe] = initial_flux/area_pipe\n",
|
||||
"\n",
|
||||
" # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n",
|
||||
" pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n",
|
||||
" p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\n",
|
||||
" v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n",
|
||||
"\n",
|
||||
" # perform the next timestep via the characteristic method\n",
|
||||
" pipe.timestep_characteristic_method()\n",
|
||||
"\n",
|
||||
" # prepare for next loop\n",
|
||||
" p_old = pipe.get_current_pressure_distribution()\n",
|
||||
" v_old = pipe.get_current_velocity_distribution()\n",
|
||||
"\n",
|
||||
" # plot some stuff\n",
|
||||
" # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n",
|
||||
" lo_00.remove()\n",
|
||||
" lo_01.remove()\n",
|
||||
" # lo_02.remove()\n",
|
||||
" # plot new pressure and velocity distribution in the pipeline\n",
|
||||
" lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,'Pa', conversion_pressure_unit),marker='.',c='blue')\n",
|
||||
" lo_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n",
|
||||
" \n",
|
||||
" fig1.suptitle(str(round(t_vec[it_pipe],2)) + '/' + str(round(t_vec[-1],2)))\n",
|
||||
" fig1.canvas.draw()\n",
|
||||
" fig1.tight_layout()\n",
|
||||
" plt.pause(0.000001)\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"fig2,axs2 = plt.subplots(2,2)\n",
|
||||
"axs2[0,0].set_title('Pressure Reservoir')\n",
|
||||
"axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit))\n",
|
||||
"axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[0,0].set_ylabel(r'$p$ [mWS]')\n",
|
||||
"axs2[0,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit))])\n",
|
||||
"\n",
|
||||
"axs2[0,1].set_title('Velocity Reservoir')\n",
|
||||
"axs2[0,1].plot(t_vec,v_boundary_res)\n",
|
||||
"axs2[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
|
||||
"axs2[0,1].set_ylim([0.9*np.min(v_boundary_res),1.1*np.max(v_boundary_res)])\n",
|
||||
"\n",
|
||||
"axs2[1,0].set_title('Pressure Turbine')\n",
|
||||
"axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit))\n",
|
||||
"axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[1,0].set_ylabel(r'$p$ [mWS]')\n",
|
||||
"axs2[1,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit))])\n",
|
||||
"\n",
|
||||
"axs2[1,1].set_title('Velocity Turbine')\n",
|
||||
"axs2[1,1].plot(t_vec,v_boundary_tur)\n",
|
||||
"axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
|
||||
"axs2[1,1].set_ylim([0.9*np.min(v_boundary_tur),1.1*np.max(v_boundary_tur)])\n",
|
||||
"\n",
|
||||
"fig2.tight_layout()\n",
|
||||
"plt.show()"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3.8.13 ('DT_Slot_3')",
|
||||
"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": "4a28055eb8a3160fa4c7e4fca69770c4e0a1add985300856aa3fcf4ce32a2c48"
|
||||
}
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 2
|
||||
}
|
||||
19
Kraftwerk/Kraftwerk_class_file.py
Normal file
19
Kraftwerk/Kraftwerk_class_file.py
Normal file
@@ -0,0 +1,19 @@
|
||||
#importing Druckrohrleitung
|
||||
import sys
|
||||
import os
|
||||
current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))
|
||||
parent = os.path.dirname(current)
|
||||
sys.path.append(parent)
|
||||
from functions.pressure_conversion import pressure_conversion
|
||||
from Turbinen.Turbinen_class_file import Francis_Turbine
|
||||
|
||||
class Kraftwerk_class:
|
||||
def __init__(self):
|
||||
self.turbines = []
|
||||
|
||||
def add_turbine(self,turbine):
|
||||
self.turbines.append(turbine)
|
||||
|
||||
def print_info(self):
|
||||
for turbine in self.turbines:
|
||||
turbine.get_info(full=True)
|
||||
105
Kraftwerk/Kraftwerk_class_test.ipynb
Normal file
105
Kraftwerk/Kraftwerk_class_test.ipynb
Normal file
@@ -0,0 +1,105 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import sys\n",
|
||||
"import os\n",
|
||||
"from Kraftwerk_class_file import Kraftwerk_class\n",
|
||||
"\n",
|
||||
"current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n",
|
||||
"parent = os.path.dirname(current)\n",
|
||||
"sys.path.append(parent)\n",
|
||||
"from functions.pressure_conversion import pressure_conversion\n",
|
||||
"from Turbinen.Turbinen_class_file import Francis_Turbine"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"[<Turbinen.Turbinen_class_file.Francis_Turbine object at 0x0000018A94FDDE80>, <Turbinen.Turbinen_class_file.Francis_Turbine object at 0x0000018A94FDDEE0>]\n",
|
||||
"Turbine has the following attributes: \n",
|
||||
"----------------------------- \n",
|
||||
"Type = Francis \n",
|
||||
"Nominal flux = 0.85 m³/s \n",
|
||||
"Nominal pressure = 108.09 mWS\n",
|
||||
"Nominal LA = 100.0 % \n",
|
||||
"Closing time = 500 s \n",
|
||||
"Current flux = -1.0 m³/s \n",
|
||||
"Current pipe pressure = -1.0 mWS \n",
|
||||
"Current LA = -1.0 % \n",
|
||||
"Simulation timestep = -1.0 s \n",
|
||||
"----------------------------- \n",
|
||||
"\n",
|
||||
"Turbine has the following attributes: \n",
|
||||
"----------------------------- \n",
|
||||
"Type = Francis \n",
|
||||
"Nominal flux = 0.85 m³/s \n",
|
||||
"Nominal pressure = 108.09 mWS\n",
|
||||
"Nominal LA = 100.0 % \n",
|
||||
"Closing time = 500 s \n",
|
||||
"Current flux = -1.0 m³/s \n",
|
||||
"Current pipe pressure = -1.0 mWS \n",
|
||||
"Current LA = -1.0 % \n",
|
||||
"Simulation timestep = -1.0 s \n",
|
||||
"----------------------------- \n",
|
||||
"\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"#Turbine\n",
|
||||
"Q_nenn = 0.85 # m³/s\n",
|
||||
"p_nenn = pressure_conversion(10.6,'bar','Pa')\n",
|
||||
"closing_time = 500 #s\n",
|
||||
"\n",
|
||||
"T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time)\n",
|
||||
"T2 = Francis_Turbine(Q_nenn,p_nenn,closing_time)\n",
|
||||
"\n",
|
||||
"KW = Kraftwerk_class()\n",
|
||||
"KW.add_turbine(T1)\n",
|
||||
"KW.add_turbine(T2)\n",
|
||||
"\n",
|
||||
"print(KW.turbines)\n",
|
||||
"\n",
|
||||
"KW.print_info()"
|
||||
]
|
||||
}
|
||||
],
|
||||
"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
|
||||
}
|
||||
@@ -2,7 +2,7 @@
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 34,
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
@@ -23,15 +23,16 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 35,
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"#define constants\n",
|
||||
"\n",
|
||||
"#Turbine\n",
|
||||
"Q_nenn = 0.85\n",
|
||||
"p_nenn,_ = pressure_conversion(10.6,'bar','Pa')\n",
|
||||
"Q_nenn = 0.85 # m³/s\n",
|
||||
"p_nenn = pressure_conversion(10.6,'bar','Pa')\n",
|
||||
"closing_time = 480. #s\n",
|
||||
"\n",
|
||||
"# physics\n",
|
||||
"g = 9.81 # gravitational acceleration [m/s²]\n",
|
||||
@@ -39,8 +40,8 @@
|
||||
"\n",
|
||||
"# define controller constants\n",
|
||||
"target_level = 8. # m\n",
|
||||
"Kp = 0.1\n",
|
||||
"Ti = 100.\n",
|
||||
"Kp = 0.01\n",
|
||||
"Ti = 3600.\n",
|
||||
"deadband_range = 0.05 # m\n",
|
||||
"\n",
|
||||
"# reservoir\n",
|
||||
@@ -59,9 +60,9 @@
|
||||
"h_fict = 100\n",
|
||||
"offset_pressure = rho*g*h_fict\n",
|
||||
"\n",
|
||||
"t_max = 1e3 #s\n",
|
||||
"nt = int(1e6) # number of simulation steps of reservoir in between timesteps of pipeline \n",
|
||||
"dt = t_max/nt\n",
|
||||
"t_max = 1e4 #s\n",
|
||||
"dt = 1e-2 # simulation timestep\n",
|
||||
"nt = int(t_max//dt) # number of simulation steps of reservoir in between timesteps of pipeline \n",
|
||||
"\n",
|
||||
"t_vec = np.arange(0,nt+1,1)*dt\n",
|
||||
"\n"
|
||||
@@ -69,25 +70,24 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 36,
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# create objects\n",
|
||||
"\n",
|
||||
"V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,dt)\n",
|
||||
"V.set_steady_state(initial_influx,initial_level,initial_pressure_unit,conversion_pressure_unit)\n",
|
||||
"V.set_steady_state(initial_influx,initial_level,conversion_pressure_unit)\n",
|
||||
"\n",
|
||||
"T1 = Francis_Turbine(Q_nenn,p_nenn)\n",
|
||||
"T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time,dt)\n",
|
||||
"T1.set_steady_state(initial_influx,p0+offset_pressure)\n",
|
||||
"T1.set_closing_time(500)\n",
|
||||
"\n",
|
||||
"Pegelregler = PI_controller_class(target_level,deadband_range,Kp,Ti,dt)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 37,
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
@@ -96,12 +96,12 @@
|
||||
"LA_soll_vec = np.full(nt+1,T1.LA)\n",
|
||||
"Q_vec = np.full(nt+1,initial_influx)\n",
|
||||
"\n",
|
||||
"Pegelregler.control_variable = T1.LA"
|
||||
"Pegelregler.control_variable = T1.get_current_LA()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 38,
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
@@ -109,106 +109,105 @@
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"0.0\n",
|
||||
"10.0\n",
|
||||
"20.0\n",
|
||||
"30.0\n",
|
||||
"40.0\n",
|
||||
"50.0\n",
|
||||
"60.0\n",
|
||||
"70.0\n",
|
||||
"80.0\n",
|
||||
"90.0\n",
|
||||
"100.0\n",
|
||||
"110.0\n",
|
||||
"120.0\n",
|
||||
"130.0\n",
|
||||
"140.0\n",
|
||||
"150.0\n",
|
||||
"160.0\n",
|
||||
"170.0\n",
|
||||
"180.0\n",
|
||||
"190.0\n",
|
||||
"200.0\n",
|
||||
"210.0\n",
|
||||
"220.0\n",
|
||||
"230.0\n",
|
||||
"240.0\n",
|
||||
"250.0\n",
|
||||
"260.0\n",
|
||||
"270.0\n",
|
||||
"280.0\n",
|
||||
"290.0\n",
|
||||
"300.0\n",
|
||||
"310.0\n",
|
||||
"320.0\n",
|
||||
"330.0\n",
|
||||
"340.0\n",
|
||||
"350.0\n",
|
||||
"360.0\n",
|
||||
"370.0\n",
|
||||
"380.0\n",
|
||||
"390.0\n",
|
||||
"400.0\n",
|
||||
"410.0\n",
|
||||
"420.0\n",
|
||||
"430.0\n",
|
||||
"440.0\n",
|
||||
"450.0\n",
|
||||
"460.0\n",
|
||||
"470.0\n",
|
||||
"480.0\n",
|
||||
"490.0\n",
|
||||
"500.0\n",
|
||||
"510.0\n",
|
||||
"520.0\n",
|
||||
"530.0\n",
|
||||
"540.0\n",
|
||||
"550.0\n",
|
||||
"560.0\n",
|
||||
"570.0\n",
|
||||
"580.0\n",
|
||||
"590.0\n",
|
||||
"600.0\n",
|
||||
"610.0\n",
|
||||
"620.0\n",
|
||||
"630.0\n",
|
||||
"640.0\n",
|
||||
"650.0\n",
|
||||
"660.0\n",
|
||||
"670.0\n",
|
||||
"680.0\n",
|
||||
"690.0\n",
|
||||
"700.0\n",
|
||||
"710.0\n",
|
||||
"720.0\n",
|
||||
"730.0\n",
|
||||
"740.0\n",
|
||||
"750.0\n",
|
||||
"760.0\n",
|
||||
"770.0\n",
|
||||
"780.0\n",
|
||||
"790.0\n",
|
||||
"800.0\n",
|
||||
"810.0\n",
|
||||
"820.0\n",
|
||||
"830.0\n",
|
||||
"840.0\n",
|
||||
"850.0\n",
|
||||
"860.0\n",
|
||||
"870.0\n",
|
||||
"880.0\n",
|
||||
"890.0\n",
|
||||
"900.0\n",
|
||||
"910.0\n",
|
||||
"920.0\n",
|
||||
"930.0\n",
|
||||
"940.0\n",
|
||||
"950.0\n",
|
||||
"960.0\n",
|
||||
"970.0\n",
|
||||
"980.0\n",
|
||||
"990.0\n",
|
||||
"1000.0\n"
|
||||
"1000.0\n",
|
||||
"1100.0\n",
|
||||
"1200.0\n",
|
||||
"1300.0\n",
|
||||
"1400.0\n",
|
||||
"1500.0\n",
|
||||
"1600.0\n",
|
||||
"1700.0\n",
|
||||
"1800.0\n",
|
||||
"1900.0\n",
|
||||
"2000.0\n",
|
||||
"2100.0\n",
|
||||
"2200.0\n",
|
||||
"2300.0\n",
|
||||
"2400.0\n",
|
||||
"2500.0\n",
|
||||
"2600.0\n",
|
||||
"2700.0\n",
|
||||
"2800.0\n",
|
||||
"2900.0\n",
|
||||
"3000.0\n",
|
||||
"3100.0\n",
|
||||
"3200.0\n",
|
||||
"3300.0\n",
|
||||
"3400.0\n",
|
||||
"3500.0\n",
|
||||
"3600.0\n",
|
||||
"3700.0\n",
|
||||
"3800.0\n",
|
||||
"3900.0\n",
|
||||
"4000.0\n",
|
||||
"4100.0\n",
|
||||
"4200.0\n",
|
||||
"4300.0\n",
|
||||
"4400.0\n",
|
||||
"4500.0\n",
|
||||
"4600.0\n",
|
||||
"4700.0\n",
|
||||
"4800.0\n",
|
||||
"4900.0\n",
|
||||
"5000.0\n",
|
||||
"5100.0\n",
|
||||
"5200.0\n",
|
||||
"5300.0\n",
|
||||
"5400.0\n",
|
||||
"5500.0\n",
|
||||
"5600.0\n",
|
||||
"5700.0\n",
|
||||
"5800.0\n",
|
||||
"5900.0\n",
|
||||
"6000.0\n",
|
||||
"6100.0\n",
|
||||
"6200.0\n",
|
||||
"6300.0\n",
|
||||
"6400.0\n",
|
||||
"6500.0\n",
|
||||
"6600.0\n",
|
||||
"6700.0\n",
|
||||
"6800.0\n",
|
||||
"6900.0\n",
|
||||
"7000.0\n",
|
||||
"7100.0\n",
|
||||
"7200.0\n",
|
||||
"7300.0\n",
|
||||
"7400.0\n",
|
||||
"7500.0\n",
|
||||
"7600.0\n",
|
||||
"7700.0\n",
|
||||
"7800.0\n",
|
||||
"7900.0\n",
|
||||
"8000.0\n",
|
||||
"8100.0\n",
|
||||
"8200.0\n",
|
||||
"8300.0\n",
|
||||
"8400.0\n",
|
||||
"8500.0\n",
|
||||
"8600.0\n",
|
||||
"8700.0\n",
|
||||
"8800.0\n",
|
||||
"8900.0\n",
|
||||
"9000.0\n",
|
||||
"9100.0\n",
|
||||
"9200.0\n",
|
||||
"9300.0\n",
|
||||
"9400.0\n",
|
||||
"9500.0\n",
|
||||
"9600.0\n",
|
||||
"9700.0\n",
|
||||
"9800.0\n",
|
||||
"9900.0\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
@@ -220,31 +219,31 @@
|
||||
" if np.mod(i,1e4) == 0:\n",
|
||||
" print(t_vec[i])\n",
|
||||
"\n",
|
||||
" if t_vec[i] == 0.4*np.max(t_vec):\n",
|
||||
" V.influx = 0\n",
|
||||
" if i == 0.4*(nt+1):\n",
|
||||
" V.set_influx(0.)\n",
|
||||
"\n",
|
||||
" p = rho*g*V.level-0.5*rho*(V.outflux_vel)**2\n",
|
||||
"\n",
|
||||
" LA_soll = Pegelregler.get_control_variable(V.level)\n",
|
||||
" T1.change_LA(LA_soll,dt)\n",
|
||||
" p = V.get_current_pressure()\n",
|
||||
" Pegelregler.update_control_variable(V.level)\n",
|
||||
" LA_soll = Pegelregler.get_current_control_variable()\n",
|
||||
" T1.update_LA(LA_soll)\n",
|
||||
" T1.set_pressure(p+offset_pressure)\n",
|
||||
" LA_soll_vec[i] = LA_soll\n",
|
||||
" LA_ist_vec[i] = T1.LA\n",
|
||||
" Q_vec[i] = T1.get_Q(p+offset_pressure)\n",
|
||||
" LA_ist_vec[i] = T1.get_current_LA()\n",
|
||||
" Q_vec[i] = T1.get_current_Q()\n",
|
||||
"\n",
|
||||
" V.pressure = p\n",
|
||||
" V.outflux_vel = 1/V.area_outflux*Q_vec[i]\n",
|
||||
" \n",
|
||||
" V.set_outflux(Q_vec[i])\n",
|
||||
"\n",
|
||||
" V.e_RK_4() \n",
|
||||
" V.level = V.update_level(V.timestep) \n",
|
||||
" V.set_volume() \n",
|
||||
" level_vec[i] = V.level \n",
|
||||
" V.timestep_reservoir_evolution() \n",
|
||||
" \n",
|
||||
" level_vec[i] = V.get_current_level()\n",
|
||||
" \n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 39,
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
@@ -257,7 +256,7 @@
|
||||
"axs1[0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs1[0].set_ylabel(r'$h$ [$\\mathrm{m}$]')\n",
|
||||
"axs1[0].plot(t_vec,level_vec)\n",
|
||||
"axs1[0].set_ylim([0.85*initial_level,1.05*initial_level])\n",
|
||||
"axs1[0].set_ylim([0*initial_level,1.5*initial_level])\n",
|
||||
"axs1[1].set_title('Flux')\n",
|
||||
"axs1[1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs1[1].set_ylabel(r'$Q$ [$\\mathrm{m} / \\mathrm{s}^3$]')\n",
|
||||
@@ -265,51 +264,33 @@
|
||||
"axs1[1].set_ylim([0,2*initial_influx])\n",
|
||||
"axs1[2].set_title('LA')\n",
|
||||
"axs1[2].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs1[2].set_ylabel(r'$LA$ [\\%]')\n",
|
||||
"axs1[2].set_ylabel(r'$LA$ [%]')\n",
|
||||
"axs1[2].plot(t_vec,LA_soll_vec)\n",
|
||||
"axs1[2].plot(t_vec,LA_ist_vec)\n",
|
||||
"axs1[2].set_ylim([0,1])\n",
|
||||
"fig1.tight_layout()\n",
|
||||
"fig1.show()\n",
|
||||
"plt.pause(1)"
|
||||
"fig1.show()\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 40,
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[<matplotlib.lines.Line2D at 0x26263b78be0>]"
|
||||
"[<matplotlib.lines.Line2D at 0x1caf15caca0>]"
|
||||
]
|
||||
},
|
||||
"execution_count": 40,
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig2 = plt.figure()\n",
|
||||
"plt.plot(t_vec,Pegelregler.error_history[1:])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 41,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"[8. 8. 8. ... 7.21126138 7.21126138 7.21126138]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"print(level_vec[:])"
|
||||
"plt.plot(t_vec,Pegelregler.get_error_history())"
|
||||
]
|
||||
}
|
||||
],
|
||||
@@ -2,13 +2,13 @@ import numpy as np
|
||||
#based on https://en.wikipedia.org/wiki/PID_controller#Discrete_implementation
|
||||
|
||||
def trap_int(vec,timestep):
|
||||
# numerical integration via the trapeziod rule to calculate the performance parameters
|
||||
l = np.size(vec)
|
||||
int = 0
|
||||
for i in range(l-1):
|
||||
int = int + (vec[i]+vec[i+1])/2*timestep
|
||||
return int
|
||||
|
||||
|
||||
def ISE_fun(error_history,timestep):
|
||||
# calcuate the integral of square error
|
||||
e = np.array(error_history)
|
||||
@@ -74,61 +74,50 @@ class P_controller_class:
|
||||
|
||||
|
||||
class PI_controller_class:
|
||||
def __init__(self,setpoint,deadband,proportionality_constant,Ti, timestep):
|
||||
# init
|
||||
def __init__(self,setpoint,deadband,proportionality_constant,Ti,timestep,lower_limit=0.,upper_limit=1.):
|
||||
self.SP = setpoint
|
||||
self.db = deadband
|
||||
self.Kp = proportionality_constant
|
||||
self.Ti = Ti
|
||||
self.Ti = Ti # integration time
|
||||
self.dt = timestep
|
||||
self.error_history = [0]
|
||||
# use a list to be able to append more easily - will get converted to np.array when needed
|
||||
self.error_history = [0]
|
||||
|
||||
self.cv_lower_limit = 0 # default
|
||||
self.cv_upper_limit = +1 # default
|
||||
self.control_variable = -99
|
||||
|
||||
self.cv_lower_limit = lower_limit # limits for the controll variable
|
||||
self.cv_upper_limit = upper_limit # limits for the controll variable
|
||||
|
||||
def set_control_variable_limits(self,lower_limit,upper_limit):
|
||||
self.cv_lower_limit = lower_limit
|
||||
self.cv_upper_limit = upper_limit
|
||||
# setter
|
||||
def set_setpoint(self,setpoint):
|
||||
self.SP = setpoint
|
||||
|
||||
def calculate_error(self,process_variable):
|
||||
self.error = process_variable-self.SP
|
||||
self.error_history.append(self.error)
|
||||
def set_control_variable(self,control_variable, display_warning=True):
|
||||
if display_warning == True and self.control_variable != -99:
|
||||
print('WARNING! You are setting the control variable of the PI controller manually \
|
||||
and are not using the .update_controll_variable() method')
|
||||
self.control_variable = control_variable
|
||||
|
||||
def get_control_variable(self,process_variable):
|
||||
|
||||
self.calculate_error(process_variable)
|
||||
|
||||
cv = self.control_variable
|
||||
Kp = self.Kp
|
||||
Ti = self.Ti
|
||||
dt = self.dt
|
||||
|
||||
e0 = self.error_history[-1]
|
||||
e1 = self.error_history[-2]
|
||||
if abs(self.error) > self.db:
|
||||
new_control = cv+Kp*(e0-e1)+dt/Ti*e0
|
||||
else:
|
||||
new_control = cv
|
||||
|
||||
if new_control < self.cv_lower_limit:
|
||||
new_control = self.cv_lower_limit
|
||||
|
||||
if new_control > self.cv_upper_limit:
|
||||
new_control = self.cv_upper_limit
|
||||
self.control_variable = new_control
|
||||
# getter
|
||||
def get_current_control_variable(self):
|
||||
return self.control_variable
|
||||
|
||||
def get_error_history(self):
|
||||
return self.error_history[1:]
|
||||
|
||||
def get_performance_indicators(self,ISE=True,IAE=True,ITSE=True,ITAE=True):
|
||||
# calculate and return the performance indicators of the error history
|
||||
ise = np.nan
|
||||
iae = np.nan
|
||||
itse = np.nan
|
||||
itae = np.nan
|
||||
|
||||
# self.error_history[1:] because the first value of the error history is set to [0]
|
||||
# to avoid special case handling in the calculation of the controll variable
|
||||
if ISE == True:
|
||||
# to avoid special case handling in the calculation of the control variable
|
||||
if ISE == True:
|
||||
ise = ISE_fun(self.error_history[1:],self.dt)
|
||||
if IAE == True:
|
||||
if IAE == True:
|
||||
iae = IAE_fun(self.error_history[1:],self.dt)
|
||||
if ITSE == True:
|
||||
itse = ITSE_fun(self.error_history[1:],self.dt)
|
||||
@@ -137,4 +126,58 @@ class PI_controller_class:
|
||||
|
||||
return ise,iae,itse,itae
|
||||
|
||||
def get_info(self):
|
||||
new_line = '\n'
|
||||
# :<10 pads the self.value to be 10 characters wide
|
||||
print_str = (f"Turbine has the following attributes: {new_line}"
|
||||
f"----------------------------- {new_line}"
|
||||
f"Type = PI Controller {new_line}"
|
||||
f"Setpoint = {self.SP:<10} {new_line}"
|
||||
f"Deadband = {self.db:<10} {new_line}"
|
||||
f"Proportionality constant = {self.Kp:<10} {new_line}"
|
||||
f"Integration time = {self.Ti:<10} [s] {new_line}"
|
||||
f"Current control variable = {round(self.control_variable,3):<10} {new_line}"
|
||||
f"Lower limit CV = {self.cv_lower_limit:<10} {new_line}"
|
||||
f"Upper limit CV = {self.cv_upper_limit:<10} {new_line}"
|
||||
f"Simulation timestep = {self.dt:<10} [s] {new_line}"
|
||||
f"----------------------------- {new_line}")
|
||||
|
||||
print(print_str)
|
||||
|
||||
# methods
|
||||
def calculate_error(self,process_variable):
|
||||
# calculate the error and expand the err history
|
||||
self.error = process_variable-self.SP
|
||||
self.error_history.append(self.error)
|
||||
|
||||
def update_control_variable(self,process_variable):
|
||||
# calculate the current control variable and make sure it does not exceed the limits
|
||||
self.calculate_error(process_variable)
|
||||
|
||||
# initialize some variables
|
||||
cv = self.control_variable
|
||||
Kp = self.Kp
|
||||
Ti = self.Ti
|
||||
dt = self.dt
|
||||
|
||||
e0 = self.error_history[-1]
|
||||
e1 = self.error_history[-2]
|
||||
|
||||
# test if the error exceeds the deadband range
|
||||
# only if that is the case, change control variable
|
||||
if abs(self.error) > self.db:
|
||||
new_control = cv+Kp*(e0-e1)+dt/Ti*e0
|
||||
else:
|
||||
new_control = cv
|
||||
|
||||
# ensure that the controll variable stays within the predefined limits
|
||||
if new_control < self.cv_lower_limit:
|
||||
new_control = self.cv_lower_limit
|
||||
if new_control > self.cv_upper_limit:
|
||||
new_control = self.cv_upper_limit
|
||||
|
||||
# set the control variable attribute
|
||||
self.set_control_variable(new_control,display_warning=False)
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -1,3 +1,4 @@
|
||||
from time import time
|
||||
import numpy as np
|
||||
#importing pressure conversion function
|
||||
import sys
|
||||
@@ -8,33 +9,126 @@ sys.path.append(parent)
|
||||
from functions.pressure_conversion import pressure_conversion
|
||||
|
||||
class Francis_Turbine:
|
||||
def __init__(self, Q_nenn,p_nenn):
|
||||
self.Q_n = Q_nenn
|
||||
self.p_n = p_nenn
|
||||
self.LA_n = 1. # 100%
|
||||
h,_ = pressure_conversion(p_nenn,'Pa','MWs')
|
||||
self.A = Q_nenn/(np.sqrt(2*9.81*h)*0.98)
|
||||
# units
|
||||
# make sure that units and print units are the same
|
||||
# units are used to label graphs and print units are used to have a bearable format when using pythons print()
|
||||
density_unit = r'$\mathrm{kg}/\mathrm{m}^3$'
|
||||
flux_unit = r'$\mathrm{m}^3/\mathrm{s}$'
|
||||
LA_unit = '%'
|
||||
pressure_unit = 'Pa'
|
||||
time_unit = 's'
|
||||
velocity_unit = r'$\mathrm{m}/\mathrm{s}$'
|
||||
volume_unit = r'$\mathrm{m}^3$'
|
||||
|
||||
def set_LA(self,LA):
|
||||
density_unit_print = 'kg/m³'
|
||||
flux_unit_print = 'm³/s'
|
||||
LA_unit_print = '%'
|
||||
pressure_unit_print = 'mWS'
|
||||
time_unit_print = 's'
|
||||
velocity_unit_print = 'm/s'
|
||||
volume_unit_print = 'm³'
|
||||
|
||||
g = 9.81 # m/s² gravitational acceleration
|
||||
|
||||
# init
|
||||
def __init__(self, Q_nenn,p_nenn,t_closing=-1.,timestep=-1.):
|
||||
self.Q_n = Q_nenn # nominal flux
|
||||
self.p_n = p_nenn # nominal pressure
|
||||
self.LA_n = 1. # 100% # nominal Leitapparatöffnung
|
||||
h = pressure_conversion(p_nenn,'Pa','MWs') # nominal pressure in terms of hydraulic head
|
||||
self.A = Q_nenn/(np.sqrt(2*self.g*h)*0.98) # Ersatzfläche
|
||||
|
||||
self.dt = timestep # simulation timestep
|
||||
self.t_c = t_closing # closing time
|
||||
self.d_LA_max_dt = 1/t_closing # maximal change of LA per second
|
||||
|
||||
# initialize for get_info() - parameters will be converted to display -1 if not overwritten
|
||||
self.p = pressure_conversion(-1,self.pressure_unit_print,self.pressure_unit)
|
||||
self.Q = -1.
|
||||
self.LA = -0.01
|
||||
|
||||
|
||||
# setter
|
||||
def set_LA(self,LA,display_warning=True):
|
||||
# set Leitapparatöffnung
|
||||
self.LA = LA
|
||||
# warn user, that the .set_LA() method should not be used ot set LA manually
|
||||
if display_warning == True:
|
||||
print('Consider using the .update_LA() method instead of setting LA manually')
|
||||
|
||||
def get_Q(self,p):
|
||||
self.Q = self.Q_n*(self.LA/self.LA_n)*np.sqrt(p/self.p_n)
|
||||
def set_timestep(self,timestep,display_warning=True):
|
||||
# set Leitapparatöffnung
|
||||
self.dt = time
|
||||
# warn user, that the .set_LA() method should not be used ot set LA manually
|
||||
if display_warning == True:
|
||||
print('WARNING: You are changing the timestep of the turbine simulation. This has implications on the simulated closing speed!')
|
||||
|
||||
def set_pressure(self,pressure):
|
||||
# set pressure in front of the turbine
|
||||
self.p = pressure
|
||||
|
||||
#getter
|
||||
def get_current_Q(self):
|
||||
# return the flux through the turbine, based on the current pressure in front
|
||||
# of the turbine and the Leitapparatöffnung
|
||||
if self.p < 0:
|
||||
self.Q = 0
|
||||
else:
|
||||
self.Q = self.Q_n*(self.LA/self.LA_n)*np.sqrt(self.p/self.p_n)
|
||||
return self.Q
|
||||
|
||||
def set_closing_time(self,t_closing):
|
||||
self.t_c = t_closing
|
||||
self.d_LA_max_dt = 1/t_closing
|
||||
def get_current_LA(self):
|
||||
return self.LA
|
||||
|
||||
def change_LA(self,LA_soll,timestep):
|
||||
LA_diff = self.LA-LA_soll
|
||||
LA_diff_max = self.d_LA_max_dt*timestep
|
||||
if abs(LA_diff) > LA_diff_max:
|
||||
LA_diff = np.sign(LA_diff)*LA_diff_max
|
||||
self.LA = self.LA-LA_diff
|
||||
def get_info(self, full = False):
|
||||
new_line = '\n'
|
||||
p = pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_print)
|
||||
p_n = pressure_conversion(self.p_n,self.pressure_unit,self.pressure_unit_print)
|
||||
|
||||
|
||||
if full == True:
|
||||
# :<10 pads the self.value to be 10 characters wide
|
||||
print_str = (f"Turbine has the following attributes: {new_line}"
|
||||
f"----------------------------- {new_line}"
|
||||
f"Type = Francis {new_line}"
|
||||
f"Nominal flux = {self.Q_n:<10} {self.flux_unit_print} {new_line}"
|
||||
f"Nominal pressure = {round(p_n,3):<10} {self.pressure_unit_print}{new_line}"
|
||||
f"Nominal LA = {self.LA_n*100:<10} {self.LA_unit_print} {new_line}"
|
||||
f"Closing time = {self.t_c:<10} {self.time_unit_print} {new_line}"
|
||||
f"Current flux = {self.Q:<10} {self.flux_unit_print} {new_line}"
|
||||
f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}"
|
||||
f"Current LA = {self.LA*100:<10} {self.LA_unit_print} {new_line}"
|
||||
f"Simulation timestep = {self.dt:<10} {self.time_unit_print} {new_line}"
|
||||
f"----------------------------- {new_line}")
|
||||
else:
|
||||
# :<10 pads the self.value to be 10 characters wide
|
||||
print_str = (f"The current attributes are: {new_line}"
|
||||
f"----------------------------- {new_line}"
|
||||
f"Current flux = {self.Q:<10} {self.flux_unit_print} {new_line}"
|
||||
f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}"
|
||||
f"Current LA = {self.LA*100:<10} {self.LA_unit_print} {new_line}"
|
||||
f"----------------------------- {new_line}")
|
||||
|
||||
print(print_str)
|
||||
|
||||
# methods
|
||||
def update_LA(self,LA_soll):
|
||||
# update the Leitappartöffnung and consider the restrictions of the closing time of the turbine
|
||||
LA_diff = self.LA-LA_soll # calculate the difference to the target LA
|
||||
LA_diff_max = self.d_LA_max_dt*self.dt # calculate the maximum change in LA based on the given timestep
|
||||
LA_diff = np.sign(LA_diff)*np.min(np.abs([LA_diff,LA_diff_max])) # calulate the correct change in LA
|
||||
|
||||
LA_new = self.LA-LA_diff
|
||||
if LA_new < 0.:
|
||||
LA_new = 0.
|
||||
elif LA_new > 1.:
|
||||
LA_new = 1.
|
||||
self.set_LA(LA_new,display_warning=False)
|
||||
|
||||
def set_steady_state(self,ss_flux,ss_pressure):
|
||||
# calculate and set steady state LA, that allows the flow of ss_flux at ss_pressure through the
|
||||
# turbine at the steady state LA
|
||||
ss_LA = self.LA_n*ss_flux/self.Q_n*np.sqrt(self.p_n/ss_pressure)
|
||||
self.set_LA(ss_LA)
|
||||
if ss_LA < 0 or ss_LA > 1:
|
||||
print('LA out of range')
|
||||
raise Exception('LA out of range [0;1]')
|
||||
self.set_LA(ss_LA,display_warning=False)
|
||||
|
||||
File diff suppressed because one or more lines are too long
266
Turbinen/Turbinen_test_Kennlinie.ipynb
Normal file
266
Turbinen/Turbinen_test_Kennlinie.ipynb
Normal file
File diff suppressed because one or more lines are too long
90
Turbinen/messy.ipynb
Normal file
90
Turbinen/messy.ipynb
Normal file
@@ -0,0 +1,90 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"from Turbinen_class_file import Francis_Turbine\n",
|
||||
"from mpl_toolkits import mplot3d\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"%matplotlib widget\n",
|
||||
"\n",
|
||||
"#importing pressure conversion function\n",
|
||||
"import sys\n",
|
||||
"import os\n",
|
||||
"current = os.path.dirname(os.path.realpath('messy.ipynb'))\n",
|
||||
"parent = os.path.dirname(current)\n",
|
||||
"sys.path.append(parent)\n",
|
||||
"from functions.pressure_conversion import pressure_conversion"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"The current attributes are: \n",
|
||||
"----------------------------- \n",
|
||||
"Current flux = -1.0 m³/s \n",
|
||||
"Current pipe pressure = -1.0 mWS \n",
|
||||
"Current LA = -1.0 % \n",
|
||||
"----------------------------- \n",
|
||||
"\n",
|
||||
"The current attributes are: \n",
|
||||
"----------------------------- \n",
|
||||
"Current flux = -1.0 m³/s \n",
|
||||
"Current pipe pressure = -1.0 mWS \n",
|
||||
"Current LA = -1.0 % \n",
|
||||
"----------------------------- \n",
|
||||
"\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"Q_nenn = 0.85\n",
|
||||
"p_nenn = pressure_conversion(10.6,'bar','Pa')\n",
|
||||
"Untertweng1 = Francis_Turbine(Q_nenn,p_nenn)\n",
|
||||
"Untertweng2 = Francis_Turbine(Q_nenn,p_nenn)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"turbines = [Untertweng1,Untertweng2]\n",
|
||||
"for turbine in turbines:\n",
|
||||
" turbine.get_info()"
|
||||
]
|
||||
}
|
||||
],
|
||||
"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
|
||||
}
|
||||
236
Untertweng.ipynb
236
Untertweng.ipynb
@@ -24,8 +24,9 @@
|
||||
"#define constants\n",
|
||||
"\n",
|
||||
"#Turbine\n",
|
||||
"Q_nenn = 0.85\n",
|
||||
"p_nenn,_ = pressure_conversion(10.6,'bar','Pa')\n",
|
||||
"Q_nenn = 0.85 # m³/s\n",
|
||||
"p_nenn = pressure_conversion(10.6,'bar','Pa')\n",
|
||||
"closing_time = 5 #s\n",
|
||||
"\n",
|
||||
"# physics\n",
|
||||
"g = 9.81 # gravitational acceleration [m/s²]\n",
|
||||
@@ -37,64 +38,40 @@
|
||||
"A_pipe = D**2/4*np.pi # pipeline area\n",
|
||||
"h_pipe = 105 # hydraulic head without reservoir [m] \n",
|
||||
"alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n",
|
||||
"n = 50 # number of pipe segments in discretization\n",
|
||||
"# consider replacing Q0 with a vector be be more flexible in initial conditions\n",
|
||||
"# Q0 = Q_nenn # initial flow in whole pipe [m³/s]\n",
|
||||
"# v0 = Q0/A_pipe # initial flow velocity [m/s]\n",
|
||||
"n = 50 # number of pipe segments in discretization # initial flow velocity [m/s]\n",
|
||||
"f_D = 0.014 # Darcy friction factor\n",
|
||||
"c = 500. # propagation velocity of the pressure wave [m/s]\n",
|
||||
"# consider prescribing a total simulation time and deducting the number of timesteps from that\n",
|
||||
"nt = 3000 # number of time steps after initial conditions\n",
|
||||
"nt = 2500 # number of time steps after initial conditions\n",
|
||||
"\n",
|
||||
"# derivatives of the pipeline constants\n",
|
||||
"dx = L/n # length of each pipe segment\n",
|
||||
"dt = dx/c # timestep according to method of characterisitics\n",
|
||||
"nn = n+1 # number of nodes\n",
|
||||
"initial_level = 8. # water level in upstream reservoir [m]\n",
|
||||
"# p0 = rho*g*initial_level-v0**2*rho/2\n",
|
||||
"pl_vec = np.arange(0,nn*dx,dx) # pl = pipe-length. position of the nodes on the pipeline\n",
|
||||
"pl_vec = np.arange(0,nn,1)*dx # pl = pipe-length. position of the nodes on the pipeline\n",
|
||||
"t_vec = np.arange(0,nt+1)*dt # time vector\n",
|
||||
"h_vec = np.arange(0,n+1)*h_pipe/n # hydraulic head of pipeline at each node \n",
|
||||
"# v_init = np.full(nn,Q0/(D**2/4*np.pi)) # initial velocity distribution in pipeline\n",
|
||||
"# p_init = (rho*g*(initial_level+h_vec)-v_init**2*rho/2)-(f_D*pl_vec/D*rho/2*v_init**2) # ref Wikipedia: Darcy Weisbach\n",
|
||||
"h_vec = np.arange(0,nn,1)*h_pipe/n # hydraulic head of pipeline at each node \n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# reservoir\n",
|
||||
"# replace influx by vector\n",
|
||||
"initial_influx = Q_nenn/1.1 # initial influx of volume to the reservoir [m³/s]\n",
|
||||
"# initial_outflux = Q0 # initial outflux of volume from the reservoir to the pipeline [m³/s]\n",
|
||||
"# initial_pipeline_pressure = p0 # Initial condition for the static pipeline pressure at the reservoir (= hydrostatic pressure - dynamic pressure) \n",
|
||||
"initial_flux = Q_nenn/1.1 # initial influx of volume to the reservoir [m³/s]\n",
|
||||
"initial_pressure_unit = 'Pa' # DO NOT CHANGE! for pressure conversion in print statements and plot labels \n",
|
||||
"conversion_pressure_unit = 'bar' # for pressure conversion in print statements and plot labels\n",
|
||||
"conversion_pressure_unit = 'bar' # for pressure conversion in print statements and plot labels\n",
|
||||
"area_base = 74. # total base are of the cuboid reservoir [m²] \n",
|
||||
"area_outflux = A_pipe # outlfux area of the reservoir, given by pipeline area [m²]\n",
|
||||
"critical_level_low = 0. # for yet-to-be-implemented warnings[m]\n",
|
||||
"critical_level_high = np.inf # for yet-to-be-implemented warnings[m]\n",
|
||||
"\n",
|
||||
"# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n",
|
||||
"nt_eRK4 = 1000 # number of simulation steps of reservoir in between timesteps of pipeline \n",
|
||||
"nt_eRK4 = 100 # number of simulation steps of reservoir in between timesteps of pipeline \n",
|
||||
"simulation_timestep = dt/nt_eRK4\n",
|
||||
"\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"#### Ideas for checks after constant definitions: \n",
|
||||
"\n",
|
||||
"- Check that the initial pressure is not negative:\n",
|
||||
" - may happen, if there is too little hydraulic head to create the initial flow conditions with the given friction\n",
|
||||
"<br>\n",
|
||||
"<br>\n",
|
||||
"- plausbility checks?\n",
|
||||
" - area > area_outflux ?\n",
|
||||
" - propable ranges for parameters?\n",
|
||||
" - angle and height/length fit together?\n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
@@ -104,21 +81,18 @@
|
||||
"# create objects\n",
|
||||
"\n",
|
||||
"V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,simulation_timestep)\n",
|
||||
"V.set_steady_state(initial_influx,initial_level,initial_pressure_unit,conversion_pressure_unit)\n",
|
||||
"V.set_steady_state(initial_flux,initial_level,conversion_pressure_unit)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n",
|
||||
"pipe.set_pressure_propagation_velocity(c)\n",
|
||||
"pipe.set_number_of_timesteps(nt)\n",
|
||||
"pipe.set_steady_state(initial_influx,V.level,pl_vec,h_vec,initial_pressure_unit,conversion_pressure_unit)\n",
|
||||
"pipe.set_steady_state(initial_flux,initial_level,area_base,pl_vec,h_vec)\n",
|
||||
"\n",
|
||||
"initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n",
|
||||
"\n",
|
||||
"T1 = Francis_Turbine(Q_nenn,p_nenn)\n",
|
||||
"T1.set_steady_state(initial_influx,pipe.p0[-1])\n",
|
||||
"T1.set_closing_time(30)\n",
|
||||
"\n",
|
||||
"# display the attributes of the created reservoir and pipeline object\n",
|
||||
"# V.get_info(full=True)\n",
|
||||
"# pipe.get_info()"
|
||||
"T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time,timestep=dt)\n",
|
||||
"T1.set_steady_state(initial_flux,initial_pressure_turbine)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -130,21 +104,20 @@
|
||||
"# initialization for timeloop\n",
|
||||
"\n",
|
||||
"# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n",
|
||||
"v_old = pipe.v0.copy()\n",
|
||||
"p_old = pipe.p0.copy()\n",
|
||||
"v_old = pipe.get_current_velocity_distribution()\n",
|
||||
"p_old = pipe.get_current_pressure_distribution()\n",
|
||||
"\n",
|
||||
"# prepare the vectors in which the temporal evolution of the boundary conditions are stored\n",
|
||||
" # keep in mind, that the velocity at the turbine and the pressure at the reservoir are set manually and\n",
|
||||
" # through the time evolution of the reservoir respectively \n",
|
||||
" # keep in mind, that the velocity at the turbine and the pressure at the reservoir follow from boundary conditions\n",
|
||||
" # reservoir level and flow through turbine\n",
|
||||
" # the pressure at the turbine and the velocity at the reservoir are calculated from the method of characteristics\n",
|
||||
"v_boundary_res = np.empty_like(t_vec)\n",
|
||||
"v_boundary_tur = np.empty_like(t_vec)\n",
|
||||
"p_boundary_res = np.empty_like(t_vec)\n",
|
||||
"p_boundary_tur = np.empty_like(t_vec)\n",
|
||||
"v_boundary_res = np.zeros_like(t_vec)\n",
|
||||
"v_boundary_tur = np.zeros_like(t_vec)\n",
|
||||
"p_boundary_res = np.zeros_like(t_vec)\n",
|
||||
"p_boundary_tur = np.zeros_like(t_vec)\n",
|
||||
"\n",
|
||||
"# prepare the vectors that store the temporal evolution of the level in the reservoir\n",
|
||||
"level_vec = np.full(nt+1,V.level) # level at the end of each pipeline timestep\n",
|
||||
"level_vec_2 = np.empty([nt_eRK4]) # level throughout each reservoir timestep-used for plotting and overwritten afterwards\n",
|
||||
"level_vec = np.full(nt+1,initial_level) # level at the end of each pipeline timestep\n",
|
||||
"\n",
|
||||
"# set the boundary conditions for the first timestep\n",
|
||||
"v_boundary_res[0] = v_old[0]\n",
|
||||
@@ -153,7 +126,8 @@
|
||||
"p_boundary_tur[0] = p_old[-1]\n",
|
||||
"\n",
|
||||
"LA_soll_vec = np.full_like(t_vec,T1.LA)\n",
|
||||
"LA_soll_vec[1500:]= 0\n",
|
||||
"LA_soll_vec[500:]= 0\n",
|
||||
"LA_ist_vec = np.full_like(t_vec,T1.LA)\n",
|
||||
"\n",
|
||||
"\n"
|
||||
]
|
||||
@@ -176,16 +150,11 @@
|
||||
"axs1[0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
|
||||
"axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
|
||||
"axs1[1].set_ylabel(r'$v$ [$\\mathrm{m} / \\mathrm{s}$]')\n",
|
||||
"lo_00, = axs1[0].plot(pl_vec,pressure_conversion(pipe.p_old,initial_pressure_unit, conversion_pressure_unit)[0],marker='.')\n",
|
||||
"lo_01, = axs1[1].plot(pl_vec,pipe.v_old,marker='.')\n",
|
||||
"lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.')\n",
|
||||
"lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n",
|
||||
"axs1[0].autoscale()\n",
|
||||
"axs1[1].autoscale()\n",
|
||||
"# displaying the reservoir level within each pipeline timestep\n",
|
||||
"# axs1[2].set_title('Level reservoir')\n",
|
||||
"# axs1[2].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"# axs1[2].set_ylabel(r'$h$ [m]')\n",
|
||||
"# lo_02, = axs1[2].plot(level_vec_2)\n",
|
||||
"# axs1[2].autoscale()\n",
|
||||
"\n",
|
||||
"fig1.tight_layout()\n",
|
||||
"fig1.show()\n",
|
||||
"plt.pause(1)\n"
|
||||
@@ -200,54 +169,54 @@
|
||||
"# loop through time steps of the pipeline\n",
|
||||
"for it_pipe in range(1,pipe.nt+1):\n",
|
||||
"\n",
|
||||
" if it_pipe == 250:\n",
|
||||
" V.set_influx(0.)\n",
|
||||
"\n",
|
||||
"# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n",
|
||||
" # set initial conditions for the reservoir time evolution calculted with e-RK4\n",
|
||||
" V.pressure = p_old[0]\n",
|
||||
" V.outflux_vel = v_old[0]\n",
|
||||
" V.set_pressure(p_old[0])\n",
|
||||
" V.set_outflux(v_old[0]*area_outflux)\n",
|
||||
" # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n",
|
||||
" for it_res in range(nt_eRK4):\n",
|
||||
" V.e_RK_4() # call e-RK4 to update outflux\n",
|
||||
" V.level = V.update_level(V.timestep) # \n",
|
||||
" V.set_volume() # update volume in reservoir\n",
|
||||
" level_vec_2[it_res] = V.level # save for plotting\n",
|
||||
" if (V.level < critical_level_low) or (V.level > critical_level_high): # make sure to never exceed critical levels\n",
|
||||
" i_max = it_pipe # for plotting only calculated values\n",
|
||||
" break \n",
|
||||
" level_vec[it_pipe] = V.level \n",
|
||||
" V.timestep_reservoir_evolution() \n",
|
||||
" level_vec[it_pipe] = V.get_current_level() \n",
|
||||
"\n",
|
||||
" # change the Leitapparatöffnung based on the target value\n",
|
||||
" T1.update_LA(LA_soll_vec[it_pipe])\n",
|
||||
" T1.set_pressure(p_old[-1])\n",
|
||||
"\n",
|
||||
" LA_ist_vec[it_pipe] = T1.LA\n",
|
||||
"\n",
|
||||
" # set boundary conditions for the next timestep of the characteristic method\n",
|
||||
" p_boundary_res[it_pipe] = rho*g*V.level-V.outflux_vel**2*rho/2\n",
|
||||
" v_boundary_res[it_pipe] = v_old[1]+1/(rho*c)*(p_boundary_res[it_pipe]-p_old[1])-f_D*dt/(2*D)*abs(v_old[1])*v_old[1] \\\n",
|
||||
" +dt*g*np.sin(alpha)\n",
|
||||
"\n",
|
||||
" T1.change_LA(LA_soll_vec[it_pipe],dt)\n",
|
||||
" v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_Q(p_old[-1])\n",
|
||||
" p_boundary_res[it_pipe] = V.get_current_pressure()\n",
|
||||
" v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_current_Q()\n",
|
||||
"\n",
|
||||
" # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n",
|
||||
" pipe.set_boundary_conditions_next_timestep(v_boundary_res[it_pipe],p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n",
|
||||
" p_boundary_tur[it_pipe] = pipe.p_boundary_tur\n",
|
||||
" pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n",
|
||||
" p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\n",
|
||||
" v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n",
|
||||
"\n",
|
||||
" # perform the next timestep via the characteristic method\n",
|
||||
" pipe.timestep_characteristic_method()\n",
|
||||
"\n",
|
||||
" # prepare for next loop\n",
|
||||
" p_old = pipe.get_current_pressure_distribution()\n",
|
||||
" v_old = pipe.get_current_velocity_distribution()\n",
|
||||
"\n",
|
||||
" # plot some stuff\n",
|
||||
" # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n",
|
||||
" lo_00.remove()\n",
|
||||
" lo_01.remove()\n",
|
||||
" # lo_02.remove()\n",
|
||||
" # plot new pressure and velocity distribution in the pipeline\n",
|
||||
" lo_00, = axs1[0].plot(pl_vec,pressure_conversion(pipe.p_old,initial_pressure_unit, conversion_pressure_unit)[0],marker='.',c='blue')\n",
|
||||
" lo_01, = axs1[1].plot(pl_vec,pipe.v_old,marker='.',c='blue')\n",
|
||||
" lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.',c='blue')\n",
|
||||
" lo_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n",
|
||||
" # lo_02, = axs1[2].plot(level_vec_2,c='blue')\n",
|
||||
" fig1.suptitle(str(round(t_vec[it_pipe],2))+ ' s / '+str(round(t_vec[-1],2)) + ' s' )\n",
|
||||
" fig1.canvas.draw()\n",
|
||||
" fig1.tight_layout()\n",
|
||||
" fig1.show()\n",
|
||||
" plt.pause(0.1) \n",
|
||||
"\n",
|
||||
" # prepare for next loop\n",
|
||||
" p_old = pipe.p_old\n",
|
||||
" v_old = pipe.v_old \n",
|
||||
" plt.pause(0.001) \n",
|
||||
"\n",
|
||||
" \n",
|
||||
" "
|
||||
@@ -257,35 +226,108 @@
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"0.0\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"print(V.get_current_influx())"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# plot time evolution of boundary pressure and velocity as well as the reservoir level\n",
|
||||
"\n",
|
||||
"fig2,axs2 = plt.subplots(3,2)\n",
|
||||
"axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,initial_pressure_unit, conversion_pressure_unit)[0])\n",
|
||||
"axs2[0,1].plot(t_vec,v_boundary_res)\n",
|
||||
"axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,initial_pressure_unit, conversion_pressure_unit)[0])\n",
|
||||
"axs2[1,1].plot(t_vec,v_boundary_tur)\n",
|
||||
"axs2[2,0].plot(t_vec,level_vec)\n",
|
||||
"axs2[0,0].set_title('Pressure reservoir')\n",
|
||||
"axs2[0,1].set_title('Velocity reservoir')\n",
|
||||
"axs2[1,0].set_title('Pressure turbine')\n",
|
||||
"axs2[1,1].set_title('Velocity turbine')\n",
|
||||
"axs2[2,0].set_title('Level reservoir')\n",
|
||||
"axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,initial_pressure_unit, conversion_pressure_unit))\n",
|
||||
"axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[0,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
|
||||
"\n",
|
||||
"axs2[0,1].set_title('Velocity reservoir')\n",
|
||||
"axs2[0,1].plot(t_vec,v_boundary_res)\n",
|
||||
"axs2[0,1].set_ylim(-2*Q_nenn,+2*Q_nenn)\n",
|
||||
"axs2[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
|
||||
"\n",
|
||||
"axs2[1,0].set_title('Pressure turbine')\n",
|
||||
"axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,initial_pressure_unit, conversion_pressure_unit))\n",
|
||||
"axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[1,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
|
||||
"\n",
|
||||
"axs2[1,1].set_title('Velocity turbine')\n",
|
||||
"axs2[1,1].plot(t_vec,v_boundary_tur)\n",
|
||||
"axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
|
||||
"\n",
|
||||
"axs2[2,0].set_title('Level reservoir')\n",
|
||||
"axs2[2,0].plot(t_vec,level_vec)\n",
|
||||
"axs2[2,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[2,0].set_ylabel(r'$h$ [m]')\n",
|
||||
"axs2[2,1].axis('off')\n",
|
||||
"\n",
|
||||
"axs2[2,1].set_title('LA')\n",
|
||||
"axs2[2,1].plot(t_vec,100*LA_soll_vec)\n",
|
||||
"axs2[2,1].plot(t_vec,100*LA_ist_vec)\n",
|
||||
"axs2[2,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[2,1].set_ylabel(r'$LA$ [%]')\n",
|
||||
"fig2.tight_layout()\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"The cuboid reservoir has the following attributes: \n",
|
||||
"----------------------------- \n",
|
||||
"Base area = 74.0 m² \n",
|
||||
"Outflux area = 0.636 m² \n",
|
||||
"Current level = 7.875725956447418 m\n",
|
||||
"Critical level low = 0.0 m \n",
|
||||
"Critical level high = inf m \n",
|
||||
"Volume in reservoir = -- m³ \n",
|
||||
"Current influx = 0.0 m³/s \n",
|
||||
"Current outflux = -0.1415386124341686 m³/s \n",
|
||||
"Current outflux vel = -0.222 m/s \n",
|
||||
"Current pipe pressure = 0.772 bar \n",
|
||||
"Simulation timestep = 0.0004052 s \n",
|
||||
"Density of liquid = 1000 kg/m³ \n",
|
||||
"----------------------------- \n",
|
||||
"\n",
|
||||
"9.22707730779877\n",
|
||||
"10.57842865915012\n",
|
||||
"11.92978001050147\n",
|
||||
"13.281131361852822\n",
|
||||
"14.632482713204173\n",
|
||||
"15.983834064555523\n",
|
||||
"17.335185415906874\n",
|
||||
"18.686536767258225\n",
|
||||
"20.037888118609576\n",
|
||||
"21.389239469960927\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"V.get_info(full=True)\n",
|
||||
"V.set_outflux(-10.)\n",
|
||||
"for i in range(10):\n",
|
||||
" V.level = V.update_level(10.)\n",
|
||||
" print(V.get_current_level())"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
|
||||
@@ -2,7 +2,7 @@
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
@@ -18,20 +18,29 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 23,
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"#define constants\n",
|
||||
"\n",
|
||||
"#Turbine\n",
|
||||
"Q_nenn = 0.85\n",
|
||||
"p_nenn,_ = pressure_conversion(10.6,'bar','Pa')\n",
|
||||
"Q_nenn = 0.85 # m³/s\n",
|
||||
"p_nenn = pressure_conversion(10.6,'bar','Pa')\n",
|
||||
"closing_time = 30. #s\n",
|
||||
"\n",
|
||||
"# physics\n",
|
||||
"g = 9.81 # gravitational acceleration [m/s²]\n",
|
||||
"rho = 1000. # density of water [kg/m³]\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# define controller constants\n",
|
||||
"target_level = 8. # m\n",
|
||||
"Kp = 0.1\n",
|
||||
"Ti = 7.\n",
|
||||
"deadband_range = 0.05 # m\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# pipeline\n",
|
||||
"L = 535.+478. # length of pipeline [m]\n",
|
||||
"D = 0.9 # pipe diameter [m]\n",
|
||||
@@ -42,111 +51,88 @@
|
||||
"f_D = 0.014 # Darcy friction factor\n",
|
||||
"c = 500. # propagation velocity of the pressure wave [m/s]\n",
|
||||
"# consider prescribing a total simulation time and deducting the number of timesteps from that\n",
|
||||
"nt = 1500 # number of time steps after initial conditions\n",
|
||||
"nt = 4500 # number of time steps after initial conditions\n",
|
||||
"\n",
|
||||
"# derivatives of the pipeline constants\n",
|
||||
"dx = L/n # length of each pipe segment\n",
|
||||
"dt = dx/c # timestep according to method of characterisitics\n",
|
||||
"nn = n+1 # number of nodes\n",
|
||||
"initial_level = 8. # water level in upstream reservoir [m]\n",
|
||||
"pl_vec = np.arange(0,nn*dx,dx) # pl = pipe-length. position of the nodes on the pipeline\n",
|
||||
"initial_level = target_level # water level in upstream reservoir [m]\n",
|
||||
"pl_vec = np.arange(0,nn,1)*dx # pl = pipe-length. position of the nodes on the pipeline\n",
|
||||
"t_vec = np.arange(0,nt+1)*dt # time vector\n",
|
||||
"h_vec = np.arange(0,n+1)*h_pipe/n # hydraulic head of pipeline at each node \n",
|
||||
"h_vec = np.arange(0,nn,1)*h_pipe/n # hydraulic head of pipeline at each node \n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# reservoir\n",
|
||||
"# replace influx by vector\n",
|
||||
"initial_influx = Q_nenn/1.1 # initial influx of volume to the reservoir [m³/s]\n",
|
||||
"initial_flux = Q_nenn/1.1 # initial influx of volume to the reservoir [m³/s]\n",
|
||||
"initial_pressure_unit = 'Pa' # DO NOT CHANGE! for pressure conversion in print statements and plot labels \n",
|
||||
"conversion_pressure_unit = 'bar' # for pressure conversion in print statements and plot labels\n",
|
||||
"conversion_pressure_unit = 'bar' # for pressure conversion in print statements and plot labels\n",
|
||||
"area_base = 74. # total base are of the cuboid reservoir [m²] \n",
|
||||
"area_outflux = A_pipe # outlfux area of the reservoir, given by pipeline area [m²]\n",
|
||||
"critical_level_low = 0. # for yet-to-be-implemented warnings[m]\n",
|
||||
"critical_level_high = np.inf # for yet-to-be-implemented warnings[m]\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# define controller constants\n",
|
||||
"target_level = initial_level # m\n",
|
||||
"Kp = 2\n",
|
||||
"Ti = 10\n",
|
||||
"deadband_range = 0.05 # m\n",
|
||||
"\n",
|
||||
"# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n",
|
||||
"nt_eRK4 = 1000 # number of simulation steps of reservoir in between timesteps of pipeline \n",
|
||||
"nt_eRK4 = 100 # number of simulation steps of reservoir in between timesteps of pipeline \n",
|
||||
"simulation_timestep = dt/nt_eRK4\n",
|
||||
"\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"#### Ideas for checks after constant definitions: \n",
|
||||
"\n",
|
||||
"- Check that the initial pressure is not negative:\n",
|
||||
" - may happen, if there is too little hydraulic head to create the initial flow conditions with the given friction\n",
|
||||
"<br>\n",
|
||||
"<br>\n",
|
||||
"- plausbility checks?\n",
|
||||
" - area > area_outflux ?\n",
|
||||
" - propable ranges for parameters?\n",
|
||||
" - angle and height/length fit together?\n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 24,
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# create objects\n",
|
||||
"\n",
|
||||
"V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,simulation_timestep)\n",
|
||||
"V.set_steady_state(initial_influx,initial_level,initial_pressure_unit,conversion_pressure_unit)\n",
|
||||
"V.set_steady_state(initial_flux,initial_level,conversion_pressure_unit)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n",
|
||||
"pipe.set_pressure_propagation_velocity(c)\n",
|
||||
"pipe.set_number_of_timesteps(nt)\n",
|
||||
"pipe.set_steady_state(initial_influx,V.level,pl_vec,h_vec,initial_pressure_unit,conversion_pressure_unit)\n",
|
||||
"pipe.set_steady_state(initial_flux,initial_level,area_base,pl_vec,h_vec)\n",
|
||||
"\n",
|
||||
"initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n",
|
||||
"\n",
|
||||
"T1 = Francis_Turbine(Q_nenn,p_nenn)\n",
|
||||
"T1.set_steady_state(initial_influx,pipe.p0[-1])\n",
|
||||
"T1.set_closing_time(5)\n",
|
||||
"T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time,timestep=dt)\n",
|
||||
"T1.set_steady_state(initial_flux,initial_pressure_turbine)\n",
|
||||
"\n",
|
||||
"T_in = Francis_Turbine(Q_nenn,p_nenn,closing_time/2,timestep=dt)\n",
|
||||
"T_in.set_steady_state(initial_flux,p_nenn)\n",
|
||||
"\n",
|
||||
"Pegelregler = PI_controller_class(target_level,deadband_range,Kp,Ti,dt)\n",
|
||||
"\n",
|
||||
"# display the attributes of the created reservoir and pipeline object\n",
|
||||
"# V.get_info(full=True)\n",
|
||||
"# pipe.get_info()"
|
||||
"Pegelregler.control_variable = T1.get_current_LA()\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 25,
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# initialization for timeloop\n",
|
||||
"\n",
|
||||
"# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n",
|
||||
"v_old = pipe.v0.copy()\n",
|
||||
"p_old = pipe.p0.copy()\n",
|
||||
"v_old = pipe.get_current_velocity_distribution()\n",
|
||||
"p_old = pipe.get_current_pressure_distribution()\n",
|
||||
"\n",
|
||||
"# prepare the vectors in which the temporal evolution of the boundary conditions are stored\n",
|
||||
" # keep in mind, that the velocity at the turbine and the pressure at the reservoir are set manually and\n",
|
||||
" # through the time evolution of the reservoir respectively \n",
|
||||
" # keep in mind, that the velocity at the turbine and the pressure at the reservoir follow from boundary conditions\n",
|
||||
" # reservoir level and flow through turbine\n",
|
||||
" # the pressure at the turbine and the velocity at the reservoir are calculated from the method of characteristics\n",
|
||||
"v_boundary_res = np.empty_like(t_vec)\n",
|
||||
"v_boundary_tur = np.empty_like(t_vec)\n",
|
||||
"p_boundary_res = np.empty_like(t_vec)\n",
|
||||
"p_boundary_tur = np.empty_like(t_vec)\n",
|
||||
"influx_vec = np.empty_like(t_vec)\n",
|
||||
"v_boundary_res = np.zeros_like(t_vec)\n",
|
||||
"v_boundary_tur = np.zeros_like(t_vec)\n",
|
||||
"p_boundary_res = np.zeros_like(t_vec)\n",
|
||||
"p_boundary_tur = np.zeros_like(t_vec)\n",
|
||||
"\n",
|
||||
"# prepare the vectors that store the temporal evolution of the level in the reservoir\n",
|
||||
"level_vec = np.full(nt+1,V.level) # level at the end of each pipeline timestep\n",
|
||||
"level_vec_2 = np.empty([nt_eRK4]) # level throughout each reservoir timestep-used for plotting and overwritten afterwards\n",
|
||||
"level_vec = np.full(nt+1,initial_level) # level at the end of each pipeline timestep\n",
|
||||
"\n",
|
||||
"# set the boundary conditions for the first timestep\n",
|
||||
"v_boundary_res[0] = v_old[0]\n",
|
||||
@@ -154,16 +140,19 @@
|
||||
"p_boundary_res[0] = p_old[0]\n",
|
||||
"p_boundary_tur[0] = p_old[-1]\n",
|
||||
"\n",
|
||||
"LA_soll_vec = np.full_like(t_vec,T1.LA)\n",
|
||||
"Pegelregler.control_variable = T1.LA\n",
|
||||
"LA_soll_vec = np.full_like(t_vec,T1.get_current_LA())\n",
|
||||
"LA_ist_vec = np.full_like(t_vec,T1.get_current_LA())\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n"
|
||||
"LA_soll_vec2 = np.full_like(t_vec,T_in.get_current_LA())\n",
|
||||
"LA_soll_vec2[500:1000] = 0.\n",
|
||||
"LA_soll_vec2[1000:1500] = 1. \n",
|
||||
"LA_soll_vec2[1500:2000] = 0.\n",
|
||||
"LA_soll_vec2[2000:2500] = 0.5 \n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 26,
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
@@ -179,16 +168,11 @@
|
||||
"axs1[0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
|
||||
"axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
|
||||
"axs1[1].set_ylabel(r'$v$ [$\\mathrm{m} / \\mathrm{s}$]')\n",
|
||||
"lo_00, = axs1[0].plot(pl_vec,pressure_conversion(pipe.p_old,initial_pressure_unit, conversion_pressure_unit)[0],marker='.')\n",
|
||||
"lo_01, = axs1[1].plot(pl_vec,pipe.v_old,marker='.')\n",
|
||||
"lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.')\n",
|
||||
"lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n",
|
||||
"axs1[0].autoscale()\n",
|
||||
"axs1[1].set_ylim([0,2])\n",
|
||||
"# displaying the reservoir level within each pipeline timestep\n",
|
||||
"# axs1[2].set_title('Level reservoir')\n",
|
||||
"# axs1[2].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"# axs1[2].set_ylabel(r'$h$ [m]')\n",
|
||||
"# lo_02, = axs1[2].plot(level_vec_2)\n",
|
||||
"# axs1[2].autoscale()\n",
|
||||
"axs1[1].autoscale()\n",
|
||||
"\n",
|
||||
"fig1.tight_layout()\n",
|
||||
"fig1.show()\n",
|
||||
"plt.pause(1)\n"
|
||||
@@ -196,64 +180,66 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 27,
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# loop through time steps of the pipeline\n",
|
||||
"for it_pipe in range(1,pipe.nt+1):\n",
|
||||
"\n",
|
||||
" if it_pipe == 150:\n",
|
||||
" V.influx = 0\n",
|
||||
" T_in.update_LA(LA_soll_vec2[it_pipe])\n",
|
||||
" T_in.set_pressure(p_nenn)\n",
|
||||
" V.set_influx(T_in.get_current_Q())\n",
|
||||
"\n",
|
||||
"# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n",
|
||||
" # set initial conditions for the reservoir time evolution calculted with e-RK4\n",
|
||||
" V.pressure = p_old[0]\n",
|
||||
" V.outflux_vel = v_old[0]\n",
|
||||
" V.set_pressure(p_old[0])\n",
|
||||
" V.set_outflux(v_old[0]*area_outflux)\n",
|
||||
" # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n",
|
||||
" for it_res in range(nt_eRK4):\n",
|
||||
" V.e_RK_4() # call e-RK4 to update outflux\n",
|
||||
" V.level = V.update_level(V.timestep) # \n",
|
||||
" V.set_volume() # update volume in reservoir\n",
|
||||
" level_vec_2[it_res] = V.level # save for plotting\n",
|
||||
" if (V.level < critical_level_low) or (V.level > critical_level_high): # make sure to never exceed critical levels\n",
|
||||
" break \n",
|
||||
" level_vec[it_pipe] = V.level \n",
|
||||
" V.timestep_reservoir_evolution() \n",
|
||||
" level_vec[it_pipe] = V.get_current_level() \n",
|
||||
"\n",
|
||||
" # get the control variable\n",
|
||||
" Pegelregler.update_control_variable(level_vec[it_pipe])\n",
|
||||
" LA_soll_vec[it_pipe] = Pegelregler.get_current_control_variable()\n",
|
||||
" \n",
|
||||
" # change the Leitapparatöffnung based on the target value\n",
|
||||
" T1.update_LA(LA_soll_vec[it_pipe])\n",
|
||||
" LA_ist_vec[it_pipe] = T1.get_current_LA()\n",
|
||||
"\n",
|
||||
" T1.set_pressure(p_old[-1])\n",
|
||||
" # set boundary conditions for the next timestep of the characteristic method\n",
|
||||
" p_boundary_res[it_pipe] = rho*g*V.level-V.outflux_vel**2*rho/2\n",
|
||||
" v_boundary_res[it_pipe] = v_old[1]+1/(rho*c)*(p_boundary_res[it_pipe]-p_old[1])-f_D*dt/(2*D)*abs(v_old[1])*v_old[1] \\\n",
|
||||
" +dt*g*np.sin(alpha)\n",
|
||||
"\n",
|
||||
" LA_soll_vec[it_pipe] = Pegelregler.get_control_variable(V.level)\n",
|
||||
" T1.change_LA(LA_soll_vec[it_pipe],dt)\n",
|
||||
" v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_Q(p_old[-1])\n",
|
||||
" p_boundary_res[it_pipe] = V.get_current_pressure()\n",
|
||||
" v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_current_Q()\n",
|
||||
"\n",
|
||||
" # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n",
|
||||
" pipe.set_boundary_conditions_next_timestep(v_boundary_res[it_pipe],p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n",
|
||||
" p_boundary_tur[it_pipe] = pipe.p_boundary_tur\n",
|
||||
" pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n",
|
||||
" p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\n",
|
||||
" v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" # perform the next timestep via the characteristic method\n",
|
||||
" pipe.timestep_characteristic_method()\n",
|
||||
"\n",
|
||||
" # prepare for next loop\n",
|
||||
" p_old = pipe.get_current_pressure_distribution()\n",
|
||||
" v_old = pipe.get_current_velocity_distribution()\n",
|
||||
"\n",
|
||||
" # plot some stuff\n",
|
||||
" # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n",
|
||||
" lo_00.remove()\n",
|
||||
" lo_01.remove()\n",
|
||||
" # lo_02.remove()\n",
|
||||
" # plot new pressure and velocity distribution in the pipeline\n",
|
||||
" lo_00, = axs1[0].plot(pl_vec,pressure_conversion(pipe.p_old,initial_pressure_unit, conversion_pressure_unit)[0],marker='.',c='blue')\n",
|
||||
" lo_01, = axs1[1].plot(pl_vec,pipe.v_old,marker='.',c='blue')\n",
|
||||
" lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.',c='blue')\n",
|
||||
" lo_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n",
|
||||
" # lo_02, = axs1[2].plot(level_vec_2,c='blue')\n",
|
||||
" fig1.suptitle(str(round(t_vec[it_pipe],2))+ ' s / '+str(round(t_vec[-1],2)) + ' s' )\n",
|
||||
" fig1.canvas.draw()\n",
|
||||
" fig1.tight_layout()\n",
|
||||
" fig1.show()\n",
|
||||
" plt.pause(0.1) \n",
|
||||
"\n",
|
||||
" # prepare for next loop\n",
|
||||
" p_old = pipe.p_old\n",
|
||||
" v_old = pipe.v_old \n",
|
||||
" plt.pause(0.001) \n",
|
||||
"\n",
|
||||
" \n",
|
||||
" "
|
||||
@@ -261,34 +247,44 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 28,
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# plot time evolution of boundary pressure and velocity as well as the reservoir level\n",
|
||||
"\n",
|
||||
"fig2,axs2 = plt.subplots(3,2)\n",
|
||||
"axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,initial_pressure_unit, conversion_pressure_unit)[0])\n",
|
||||
"axs2[0,1].plot(t_vec,v_boundary_res)\n",
|
||||
"axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,initial_pressure_unit, conversion_pressure_unit)[0])\n",
|
||||
"axs2[1,1].plot(t_vec,v_boundary_tur)\n",
|
||||
"axs2[2,0].plot(t_vec,level_vec)\n",
|
||||
"axs2[0,0].set_title('Pressure reservoir')\n",
|
||||
"axs2[0,1].set_title('Velocity reservoir')\n",
|
||||
"axs2[1,0].set_title('Pressure turbine')\n",
|
||||
"axs2[1,1].set_title('Velocity turbine')\n",
|
||||
"axs2[2,0].set_title('Level reservoir')\n",
|
||||
"axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,initial_pressure_unit, conversion_pressure_unit))\n",
|
||||
"axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[0,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
|
||||
"\n",
|
||||
"axs2[0,1].set_title('Velocity reservoir')\n",
|
||||
"axs2[0,1].plot(t_vec,v_boundary_res)\n",
|
||||
"axs2[0,1].set_ylim(-2*Q_nenn,+2*Q_nenn)\n",
|
||||
"axs2[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
|
||||
"\n",
|
||||
"axs2[1,0].set_title('Pressure turbine')\n",
|
||||
"axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,initial_pressure_unit, conversion_pressure_unit))\n",
|
||||
"axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[1,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
|
||||
"\n",
|
||||
"axs2[1,1].set_title('Velocity turbine')\n",
|
||||
"axs2[1,1].plot(t_vec,v_boundary_tur)\n",
|
||||
"axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
|
||||
"\n",
|
||||
"axs2[2,0].set_title('Level reservoir')\n",
|
||||
"axs2[2,0].plot(t_vec,level_vec)\n",
|
||||
"axs2[2,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[2,0].set_ylabel(r'$h$ [m]')\n",
|
||||
"axs2[2,1].axis('off')\n",
|
||||
"\n",
|
||||
"axs2[2,1].set_title('LA')\n",
|
||||
"axs2[2,1].plot(t_vec,100*LA_soll_vec)\n",
|
||||
"axs2[2,1].plot(t_vec,100*LA_ist_vec)\n",
|
||||
"axs2[2,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
|
||||
"axs2[2,1].set_ylabel(r'$LA$ [%]')\n",
|
||||
"fig2.tight_layout()\n",
|
||||
"plt.show()"
|
||||
]
|
||||
|
||||
@@ -32,7 +32,7 @@ def pa_to_atm(p):
|
||||
def pa_to_psi(p):
|
||||
return p/6894.8
|
||||
|
||||
def pressure_conversion(pressure, input_unit = 'bar', target_unit = 'Pa'):
|
||||
def pressure_conversion(pressure, input_unit = 'bar', target_unit = 'Pa', return_unit = False):
|
||||
p = pressure
|
||||
if input_unit.lower() == 'bar':
|
||||
p_pa = bar_to_pa(p)
|
||||
@@ -50,20 +50,27 @@ def pressure_conversion(pressure, input_unit = 'bar', target_unit = 'Pa'):
|
||||
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
|
||||
return_vec = [pa_to_bar(p_pa), target_unit]
|
||||
elif target_unit.lower() == 'mws':
|
||||
return pa_to_mWS(p_pa), target_unit
|
||||
return_vec = [pa_to_mWS(p_pa), target_unit]
|
||||
elif target_unit.lower() == 'torr':
|
||||
return pa_to_torr(p_pa), target_unit
|
||||
return_vec = [pa_to_torr(p_pa), target_unit]
|
||||
elif target_unit.lower() == 'atm':
|
||||
return pa_to_atm(p_pa), target_unit
|
||||
return_vec = [pa_to_atm(p_pa), target_unit]
|
||||
elif target_unit.lower() =='psi':
|
||||
return pa_to_psi(p_pa), target_unit
|
||||
return_vec = [pa_to_psi(p_pa), target_unit]
|
||||
elif target_unit.lower() == 'pa':
|
||||
return p_pa, target_unit
|
||||
return_vec = [p_pa, target_unit]
|
||||
else:
|
||||
raise Exception('Given target unit not recognised. \n Known units are: Pa, bar, mWs, Torr, atm, psi')
|
||||
|
||||
if return_unit == True:
|
||||
# return with pressure unit
|
||||
return return_vec
|
||||
else:
|
||||
# return without pressure unit
|
||||
return return_vec[0]
|
||||
|
||||
# testing_pressure_conversion
|
||||
if __name__ == '__main__':
|
||||
p = 1
|
||||
@@ -72,6 +79,6 @@ if __name__ == '__main__':
|
||||
|
||||
for input_unit in unit_dict:
|
||||
for target_unit in unit_dict:
|
||||
converted_p = pressure_conversion(p,input_unit,target_unit)
|
||||
converted_p = pressure_conversion(p,input_unit,target_unit,return_unit=False)
|
||||
print(input_unit,target_unit)
|
||||
print(converted_p)
|
||||
Reference in New Issue
Block a user