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

Merge branch 'Dev'

Browse Source
This commit is contained in:
Brantegger Georg
2022-07-25 10:28:37 +02:00
parent 84ccadc66d 9204729d0b
commit 9c70ce3d8d
16 changed files with 2270 additions and 89 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
Ausgleichsbecken/Ausgleichsbecken_class_file.py
View File

@@ -1,5 +1,4 @@
import numpy as np import numpy as np
# from Ausgleichsbecken_functions import FODE_function, get_h_halfstep, get_p_halfstep
#importing pressure conversion function #importing pressure conversion function
import sys import sys
@@ -34,6 +33,9 @@ class Ausgleichsbecken_class:
velocity_unit_print = 'm/s' velocity_unit_print = 'm/s'
volume_unit_print = 'm³' volume_unit_print = 'm³'
g = 9.81 # m/s²
rho = 1000 # kg/m³
# init # 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):
self.area = area # base area of the rectangular structure self.area = area # base area of the rectangular structure
@@ -58,16 +60,26 @@ class Ausgleichsbecken_class:
self.set_volume() self.set_volume()
def set_influx(self,influx): def set_influx(self,influx):
self.influx = influx self.influx = influx
def set_outflux(self,outflux): def set_outflux(self,outflux):
self.outflux = outflux self.outflux = outflux
self.outflux_vel = outflux/self.area_outflux self.outflux_vel = outflux/self.area_outflux
def set_pressure(self,pressure,pressure_unit,display_pressure_unit): def set_pressure(self,pressure,pressure_unit,display_pressure_unit):
self.pressure = pressure self.pressure = pressure
self.pressure_unit = pressure_unit self.pressure_unit = pressure_unit
self.pressure_unit_print = display_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
self.set_initial_level(ss_level)
self.set_influx(ss_influx)
self.set_outflux(ss_outflux)
self.set_pressure(ss_pressure,pressure_unit,display_pressure_unit)
# getter # getter
def get_info(self, full = False): def get_info(self, full = False):
new_line = '\n' new_line = '\n'
@@ -86,7 +98,7 @@ class Ausgleichsbecken_class:
f"Volume in reservoir = {self.volume:<10} {self.volume_unit_print} {new_line}" 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 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 = {self.outflux:<10} {self.flux_unit_print} {new_line}"
f"Current outflux vel = {self.outflux_vel:<10} {self.velocity_unit_print} {new_line}" f"Current outflux vel = {round(self.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"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"Simulation timestep = {self.timestep:<10} {self.time_unit_print} {new_line}"
f"----------------------------- {new_line}") f"----------------------------- {new_line}")
@@ -98,7 +110,7 @@ class Ausgleichsbecken_class:
f"Volume in reservoir = {self.volume:<10} {self.volume_unit_print} {new_line}" 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 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 = {self.outflux:<10} {self.flux_unit_print} {new_line}"
f"Current outflux vel = {self.outflux_vel:<10} {self.velocity_unit_print} {new_line}" f"Current outflux vel = {round(self.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"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}"
f"----------------------------- {new_line}") f"----------------------------- {new_line}")

60
Ausgleichsbecken/Main_Program.ipynb → Ausgleichsbecken/Ausgleichsbecken_test.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 16, "execution_count": 12,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -21,16 +21,16 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 17, "execution_count": 13,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# define constants\n", "# define constants\n",
"initial_level = 5. # m\n", "initial_level = 10. # m\n",
"initial_influx = 0.5 # m³/s\n", "initial_influx = 5. # m³/s\n",
"initial_outflux = 0. # m³/s\n", "initial_outflux = 1. # m³/s\n",
"initial_pipeline_pressure = 1\n", "initial_pipeline_pressure = 10.\n",
"initial_pressure_unit = 'bar'\n", "initial_pressure_unit = 'mWS'\n",
"conversion_pressure_unit = 'mWS'\n", "conversion_pressure_unit = 'mWS'\n",
"\n", "\n",
"area_base = 1. # m²\n", "area_base = 1. # m²\n",
@@ -41,33 +41,33 @@
"\n", "\n",
"# for while loop\n", "# for while loop\n",
"total_min_level = 0.01 # m\n", "total_min_level = 0.01 # m\n",
"total_max_time = 300 # s" "total_max_time = 1000 # s"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 18, "execution_count": 14,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"%matplotlib qt\n", "%matplotlib qt\n",
"\n", "\n",
"V = Ausgleichsbecken_class(area_base, area_outflux, critical_level_low, critical_level_high,simulation_timestep)\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_initial_level(initial_level) \n",
"V.set_influx(initial_influx)\n", "# V.set_influx(initial_influx)\n",
"V.set_outflux(initial_outflux)\n", "# V.set_outflux(initial_outflux)\n",
"\n", "# converted_pressure,_ = pressure_conversion(initial_pipeline_pressure,input_unit = initial_pressure_unit, target_unit = 'Pa')\n",
"converted_pressure, V.pressure_unit = pressure_conversion(initial_pipeline_pressure,input_unit = initial_pressure_unit, target_unit = conversion_pressure_unit)\n", "# V.pressure = converted_pressure\n",
"V.pressure = converted_pressure\n", "V.set_steady_state(initial_influx,initial_level,initial_pressure_unit,conversion_pressure_unit)\n",
"\n", "\n",
"time_vec = np.arange(0,total_max_time,simulation_timestep)\n", "time_vec = np.arange(0,total_max_time,simulation_timestep)\n",
"outflux_vec = np.empty_like(time_vec)\n", "outflux_vec = np.empty_like(time_vec)\n",
"outflux_vec[0] = initial_outflux\n", "outflux_vec[0] = V.outflux\n",
"level_vec = np.empty_like(time_vec)\n", "level_vec = np.empty_like(time_vec)\n",
"level_vec[0] = initial_level\n", "level_vec[0] = V.level\n",
"\n", "\n",
"pressure_vec = np.full_like(time_vec,converted_pressure)*((np.sin(time_vec/5)+1)*np.exp(-time_vec/50))\n", "# pressure_vec = np.full_like(time_vec,converted_pressure)*((np.sin(time_vec)+1)*np.exp(-time_vec/50))\n",
" \n", "pressure_vec = np.full_like(time_vec,V.pressure)\n",
" \n", " \n",
"i_max = -1\n", "i_max = -1\n",
"\n", "\n",
@@ -82,7 +82,15 @@
" if V.level < total_min_level:\n", " if V.level < total_min_level:\n",
" i_max = i\n", " i_max = i\n",
" break\n", " break\n",
"\n", "\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"\n", "\n",
"fig1, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1)\n", "fig1, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1)\n",
"fig1.set_figheight(10)\n", "fig1.set_figheight(10)\n",
@@ -98,16 +106,16 @@
"ax2.set_xlabel(r'$t$ ['+V.time_unit+']')\n", "ax2.set_xlabel(r'$t$ ['+V.time_unit+']')\n",
"ax2.legend()\n", "ax2.legend()\n",
"\n", "\n",
"ax3.plot(time_vec[:i_max],pressure_vec[:i_max], label='Pipeline pressure at reservoir')\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}$ ['+V.pressure_unit+']')\n", "ax3.set_ylabel(r'$p_{pipeline}$ ['+conversion_pressure_unit+']')\n",
"ax3.set_xlabel(r'$t$ ['+V.time_unit+']')\n", "ax3.set_xlabel(r'$t$ ['+V.time_unit+']')\n",
"ax3.legend()\n", "ax3.legend()\n",
"\n", "\n",
"# plt.subplots_adjust(left=0.2, bottom=0.2)\n", "# plt.subplots_adjust(left=0.2, bottom=0.2)\n",
"ax4.set_axis_off()\n", "ax4.set_axis_off()\n",
"cell_text = np.array([[initial_level, V.level_unit], \\\n", "cell_text = np.array([[level_vec[0], V.level_unit], \\\n",
" [initial_influx, V.flux_unit], \\\n", " [initial_influx, V.flux_unit], \\\n",
" [initial_outflux, V.flux_unit], \\\n", " [outflux_vec[0], V.flux_unit], \\\n",
" [simulation_timestep, V.time_unit], \\\n", " [simulation_timestep, V.time_unit], \\\n",
" [area_base, V.area_unit], \\\n", " [area_base, V.area_unit], \\\n",
" [area_outflux, V.area_unit]])\n", " [area_outflux, V.area_unit]])\n",
@@ -133,7 +141,7 @@
], ],
"metadata": { "metadata": {
"kernelspec": { "kernelspec": {
"display_name": "Python 3.8.13 ('Georg_DT_Slot3')", "display_name": "Python 3.8.13 ('DT_Slot_3')",
"language": "python", "language": "python",
"name": "python3" "name": "python3"
}, },
@@ -152,7 +160,7 @@
"orig_nbformat": 4, "orig_nbformat": 4,
"vscode": { "vscode": {
"interpreter": { "interpreter": {
"hash": "84fb123bdc47ab647d3782661abcbe80fbb79236dd2f8adf4cef30e8755eb2cd" "hash": "4a28055eb8a3160fa4c7e4fca69770c4e0a1add985300856aa3fcf4ce32a2c48"
} }
} }
}, },

16
Druckrohrleitung/Druckrohrleitung_class_file.py
View File

@@ -41,7 +41,7 @@ class Druckrohrleitung_class:
self.n_seg = number_segments self.n_seg = number_segments
self.angle = pipeline_angle self.angle = pipeline_angle
self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient
self.density = rho self.rho = rho
self.g = g self.g = g
self.dx = total_length/number_segments self.dx = total_length/number_segments
@@ -89,8 +89,8 @@ class Druckrohrleitung_class:
self.v_old = self.v0.copy() self.v_old = self.v0.copy()
self.v = np.empty_like(self.v_old) self.v = np.empty_like(self.v_old)
def set_boundary_conditions_next_timestep(self,v_reservoir,p_reservoir,v_turbine,input_unit_pressure = 'Pa'): def set_boundary_conditions_next_timestep(self,v_reservoir,p_reservoir,v_turbine):
rho = self.density rho = self.rho
c = self.c c = self.c
f_D = self.f_D f_D = self.f_D
dt = self.dt dt = self.dt
@@ -108,6 +108,14 @@ class Druckrohrleitung_class:
self.p[0] = self.p_boundary_res.copy() self.p[0] = self.p_boundary_res.copy()
self.p[-1] = self.p_boundary_tur.copy() self.p[-1] = self.p_boundary_tur.copy()
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)
self.set_initial_flow_velocity(ss_v0)
self.set_initial_pressure(ss_pressure,pressure_unit,display_pressure_unit)
# getter # getter
def get_info(self): def get_info(self):
new_line = '\n' new_line = '\n'
@@ -150,7 +158,7 @@ class Druckrohrleitung_class:
def timestep_characteristic_method(self): def timestep_characteristic_method(self):
#number of nodes #number of nodes
nn = self.n_seg+1 nn = self.n_seg+1
rho = self.density rho = self.rho
c = self.c c = self.c
f_D = self.f_D f_D = self.f_D
dt = self.dt dt = self.dt

0
Druckrohrleitung/Main_Programm.ipynb → Druckrohrleitung/Druckrohrleitung_test.ipynb
View File

140
Regler/Regler_class_file.py Normal file
View File

@@ -0,0 +1,140 @@
import numpy as np
#based on https://en.wikipedia.org/wiki/PID_controller#Discrete_implementation
def trap_int(vec,timestep):
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)
dt = timestep
ise = trap_int(e**2,dt)
return ise
def IAE_fun(error_history,timestep):
# calcuate the integral of absolute error
e = np.array(error_history)
dt = timestep
iae = trap_int(np.abs(e),dt)
return iae
def ITSE_fun(error_history,timestep):
# calcuate the integral of time multiply square error
e = np.array(error_history)
dt = timestep
n = np.size(e)
t = np.arange(0,n)*dt
itse = trap_int(t*e**2,dt)
return itse
def ITAE_fun(error_history,timestep):
# calcuate the integral of time multiply absolute error
e = np.array(error_history)
dt = timestep
n = np.size(e)
t = np.arange(0,n)*dt
itae = trap_int(np.abs(e),dt)
return itae
class P_controller_class:
# def __init__(self,setpoint,proportionality_constant):
# self.SP = setpoint
# self.Kp = proportionality_constant
# self.error_history = []
# self.control_variable = 0.1
# self.lower_limit = -0.1 # default
# self.upper_limit = +0.1 # default
# def set_control_variable_limits(self,lower_limit,upper_limit):
# self.lower_limit = lower_limit
# self.upper_limit = upper_limit
# def calculate_error(self,process_variable):
# self.error = self.SP-process_variable
# self.error_history.append(self.error)
# def get_control_variable(self):
# new_control = self.control_variable+self.Kp*(self.error_history[-1]-self.error_history[-2])
# if new_control < self.lower_limit:
# new_control = self.lower_limit
# if new_control > self.upper_limit:
# new_control = self.upper_limit
# self.control_variable = new_control
# # print(new_control)
# return new_control
def __init__(self):
pass
class PI_controller_class:
def __init__(self,setpoint,deadband,proportionality_constant,Ti, timestep):
self.SP = setpoint
self.db = deadband
self.Kp = proportionality_constant
self.Ti = Ti
self.dt = timestep
self.error_history = [0]
self.cv_lower_limit = 0 # default
self.cv_upper_limit = +1 # default
def set_control_variable_limits(self,lower_limit,upper_limit):
self.cv_lower_limit = lower_limit
self.cv_upper_limit = upper_limit
def calculate_error(self,process_variable):
self.error = process_variable-self.SP
self.error_history.append(self.error)
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
return self.control_variable
def get_performance_indicators(self,ISE=True,IAE=True,ITSE=True,ITAE=True):
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:
ise = ISE_fun(self.error_history[1:],self.dt)
if IAE == True:
iae = IAE_fun(self.error_history[1:],self.dt)
if ITSE == True:
itse = ITSE_fun(self.error_history[1:],self.dt)
if ITAE == True:
itae = ITAE_fun(self.error_history[1:],self.dt)
return ise,iae,itse,itae

343
Regler/regler_test.ipynb Normal file
View File

@@ -0,0 +1,343 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from Regler_class_file import PI_controller_class\n",
"\n",
"#importing Druckrohrleitung\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\n",
"from Turbinen.Turbinen_class_file import Francis_Turbine"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"#define constants\n",
"\n",
"#Turbine\n",
"Q_nenn = 0.85\n",
"p_nenn,_ = pressure_conversion(10.6,'bar','Pa')\n",
"\n",
"# physics\n",
"g = 9.81 # gravitational acceleration [m/s²]\n",
"rho = 1000. # density of water [kg/m³]\n",
"\n",
"# define controller constants\n",
"target_level = 8. # m\n",
"Kp = 0.1\n",
"Ti = 100.\n",
"deadband_range = 0.05 # m\n",
"\n",
"# reservoir\n",
"initial_level = target_level\n",
"initial_influx = Q_nenn/2 # 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",
"area_base = 74. # total base are of the cuboid reservoir [m²] \n",
"area_outflux = 1. # outflux 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",
"p0 = rho*g*initial_level-0.5*rho*(initial_influx/area_outflux)**2\n",
"\n",
"# offset the pressure in front of the turbine to get realisitc fluxes\n",
"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",
"\n",
"t_vec = np.arange(0,nt+1,1)*dt\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 36,
"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",
"\n",
"T1 = Francis_Turbine(Q_nenn,p_nenn)\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,
"metadata": {},
"outputs": [],
"source": [
"level_vec = np.full(nt+1,V.level)\n",
"LA_ist_vec = np.full(nt+1,T1.LA)\n",
"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"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"name": "stdout",
"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"
]
}
],
"source": [
"# time loop\n",
"\n",
"for i in range(nt+1):\n",
"\n",
" 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",
"\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",
" LA_soll_vec[i] = LA_soll\n",
" LA_ist_vec[i] = T1.LA\n",
" Q_vec[i] = T1.get_Q(p+offset_pressure)\n",
"\n",
" V.pressure = p\n",
" V.outflux_vel = 1/V.area_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",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib qt5\n",
"# time loop\n",
"\n",
"# create a figure and subplots to display the velocity and pressure distribution across the pipeline in each pipeline step\n",
"fig1,axs1 = plt.subplots(3,1)\n",
"axs1[0].set_title('Level')\n",
"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[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",
"axs1[1].plot(t_vec,Q_vec)\n",
"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].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)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x26263b78be0>]"
]
},
"execution_count": 40,
"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[:])"
]
}
],
"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
}

386
Regler/regler_test_optimierung.ipynb Normal file
View File

File diff suppressed because one or more lines are too long

40
Turbinen/Turbinen_class_file.py Normal file
View File

@@ -0,0 +1,40 @@
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 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)
def set_LA(self,LA):
self.LA = LA
def get_Q(self,p):
self.Q = self.Q_n*(self.LA/self.LA_n)*np.sqrt(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 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 set_steady_state(self,ss_flux,ss_pressure):
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')

239
Turbinen/Turbinen_test.ipynb Normal file
View File

File diff suppressed because one or more lines are too long

22
Turbinen/old/Durchflusskennlinie.csv Normal file
View File

@@ -0,0 +1,22 @@
,11.4,11.2,11,10.8,10.6,10.4,10.2,10,9.8
0,0,0,0,0,0,0,0,0,0
0.05,44.6719225,43.934144,43.3914212,43.005945,42.7411852,42.5620659,42.4351104,42.3285595,42.2124611
0.1,93.5257218,92.1813802,91.0120507,89.9819869,89.0566946,88.2030946,87.3896575,86.5865116,85.7655241
0.15,142.455373,140.502298,138.703994,137.026824,135.438371,133.907593,132.404945,130.902474,129.373898
0.2,191.35358,188.792245,186.365298,184.041241,181.789769,179.581903,177.390108,175.188376,172.952294
0.25,240.112708,236.946245,233.893698,230.92573,228.014163,225.132101,222.254034,219.355912,216.415204
0.3,288.625576,284.85976,281.187353,277.581187,274.01522,270.464644,266.905977,263.31713,259.677456
0.35,336.786234,332.429439,328.145567,323.909615,319.697669,315.487006,311.256165,306.985012,302.654777
0.4,384.490739,379.553866,374.669505,369.814802,364.967956,360.108307,355.216403,350.274048,345.264331
0.45,431.637894,426.134271,420.662881,415.202987,409.734875,404.239922,398.700655,393.100789,387.425251
0.5,478.129951,472.075209,466.032607,459.983487,453.910176,447.796055,441.625591,435.384378,429.059145
0.55,523.873268,517.285198,510.689413,504.069281,497.409128,490.694283,483.911113,477.047044,470.090565
0.6,568.778912,561.677293,554.548395,547.377555,540.151033,532.856054,525.480827,518.014558,510.447451
0.65,612.763186,605.169605,597.529525,589.830179,582.059697,574.207132,566.262474,558.216649,550.061519
0.7,655.7481,647.685753,639.558081,631.354134,623.063835,614.677994,606.188309,597.587364,588.868614
0.75,697.661758,689.155243,680.565018,671.881864,663.097416,654.204159,645.195426,636.065384,626.809013
0.8,738.438667,729.51377,720.487263,711.35157,702.099947,692.726469,683.226022,673.594278,663.827671
0.85,778.019972,768.703447,759.267942,749.707427,740.016685,730.191293,720.227602,710.122707,699.874419
0.9,816.35361,806.672962,796.856534,786.899741,776.798797,766.550685,756.153132,745.604572,734.904109
0.95,853.394385,843.377654,833.208949,822.885029,812.403437,801.762466,790.961126,779.999101,768.876705
1,889.103974,878.779525,868.287549,857.626044,846.793778,835.790258,824.615682,813.270891,801.757325
1 11.4 11.2 11 10.8 10.6 10.4 10.2 10 9.8
2 0 0 0 0 0 0 0 0 0 0
3 0.05 44.6719225 43.934144 43.3914212 43.005945 42.7411852 42.5620659 42.4351104 42.3285595 42.2124611
4 0.1 93.5257218 92.1813802 91.0120507 89.9819869 89.0566946 88.2030946 87.3896575 86.5865116 85.7655241
5 0.15 142.455373 140.502298 138.703994 137.026824 135.438371 133.907593 132.404945 130.902474 129.373898
6 0.2 191.35358 188.792245 186.365298 184.041241 181.789769 179.581903 177.390108 175.188376 172.952294
7 0.25 240.112708 236.946245 233.893698 230.92573 228.014163 225.132101 222.254034 219.355912 216.415204
8 0.3 288.625576 284.85976 281.187353 277.581187 274.01522 270.464644 266.905977 263.31713 259.677456
9 0.35 336.786234 332.429439 328.145567 323.909615 319.697669 315.487006 311.256165 306.985012 302.654777
10 0.4 384.490739 379.553866 374.669505 369.814802 364.967956 360.108307 355.216403 350.274048 345.264331
11 0.45 431.637894 426.134271 420.662881 415.202987 409.734875 404.239922 398.700655 393.100789 387.425251
12 0.5 478.129951 472.075209 466.032607 459.983487 453.910176 447.796055 441.625591 435.384378 429.059145
13 0.55 523.873268 517.285198 510.689413 504.069281 497.409128 490.694283 483.911113 477.047044 470.090565
14 0.6 568.778912 561.677293 554.548395 547.377555 540.151033 532.856054 525.480827 518.014558 510.447451
15 0.65 612.763186 605.169605 597.529525 589.830179 582.059697 574.207132 566.262474 558.216649 550.061519
16 0.7 655.7481 647.685753 639.558081 631.354134 623.063835 614.677994 606.188309 597.587364 588.868614
17 0.75 697.661758 689.155243 680.565018 671.881864 663.097416 654.204159 645.195426 636.065384 626.809013
18 0.8 738.438667 729.51377 720.487263 711.35157 702.099947 692.726469 683.226022 673.594278 663.827671
19 0.85 778.019972 768.703447 759.267942 749.707427 740.016685 730.191293 720.227602 710.122707 699.874419
20 0.9 816.35361 806.672962 796.856534 786.899741 776.798797 766.550685 756.153132 745.604572 734.904109
21 0.95 853.394385 843.377654 833.208949 822.885029 812.403437 801.762466 790.961126 779.999101 768.876705
22 1 889.103974 878.779525 868.287549 857.626044 846.793778 835.790258 824.615682 813.270891 801.757325

33
Turbinen/old/Turbinen_class_file.py Normal file
View File

@@ -0,0 +1,33 @@
from matplotlib.pyplot import fill
import numpy as np
from scipy.interpolate import interp2d
#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 Francis_turbine_class:
def __init__(self,CSV_name='Durchflusskennlinie.csv'):
csv = np.genfromtxt(CSV_name,delimiter=',')
n_rows,_ = np.shape(csv)
self.raw_csv = np.append(csv,np.zeros([n_rows,1]),axis = 1)
def extract_csv(self,CSV_pressure_unit='bar'):
ps_vec,_ = pressure_conversion(self.raw_csv[0,1:],CSV_pressure_unit,'Pa')
self.raw_ps_vec = np.flip(ps_vec)
self.raw_LA_vec = self.raw_csv[1:,0]
self.raw_Qs_mat = np.fliplr(self.raw_csv[1:,1:])/1000. # convert from l/s to m³/s
def get_Q_fun(self):
Q_fun = interp2d(self.raw_ps_vec,self.raw_LA_vec,self.raw_Qs_mat,bounds_error=False,fill_value=None)
return Q_fun

319
Turbinen/old/Turbinen_test.ipynb Normal file
View File

File diff suppressed because one or more lines are too long

278
Turbinen/old/messy.ipynb Normal file
View File

File diff suppressed because one or more lines are too long

115
Main_Programm.ipynb → Untertweng.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 11, "execution_count": 1,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -11,57 +11,62 @@
"\n", "\n",
"from functions.pressure_conversion import pressure_conversion\n", "from functions.pressure_conversion import pressure_conversion\n",
"from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class\n", "from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class\n",
"from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class" "from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class\n",
"from Turbinen.Turbinen_class_file import Francis_Turbine"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 12, "execution_count": 2,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"#define constants\n", "#define constants\n",
"\n", "\n",
"#Turbine\n",
"Q_nenn = 0.85\n",
"p_nenn,_ = pressure_conversion(10.6,'bar','Pa')\n",
"\n",
"# physics\n", "# physics\n",
"g = 9.81 # gravitational acceleration [m/s²]\n", "g = 9.81 # gravitational acceleration [m/s²]\n",
"rho = 1000. # density of water [kg/m³]\n", "rho = 1000. # density of water [kg/m³]\n",
"\n", "\n",
"# pipeline\n", "# pipeline\n",
"L = 1000. # length of pipeline [m]\n", "L = 535.+478. # length of pipeline [m]\n",
"D = 1. # pipe diameter [m]\n", "D = 0.9 # pipe diameter [m]\n",
"A_pipe = D**2/4*np.pi # pipeline area\n", "A_pipe = D**2/4*np.pi # pipeline area\n",
"h_pipe = 200 # hydraulic head without reservoir [m] \n", "h_pipe = 105 # hydraulic head without reservoir [m] \n",
"alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n", "alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n",
"n = 100 # number of pipe segments in discretization\n", "n = 50 # number of pipe segments in discretization\n",
"# consider replacing Q0 with a vector be be more flexible in initial conditions\n", "# consider replacing Q0 with a vector be be more flexible in initial conditions\n",
"Q0 = 2. # initial flow in whole pipe [m³/s]\n", "# Q0 = Q_nenn # initial flow in whole pipe [m³/s]\n",
"v0 = Q0/A_pipe # initial flow velocity [m/s]\n", "# v0 = Q0/A_pipe # initial flow velocity [m/s]\n",
"f_D = 0.1 # Darcy friction factor\n", "f_D = 0.014 # Darcy friction factor\n",
"c = 400. # propagation velocity of the pressure wave [m/s]\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", "# consider prescribing a total simulation time and deducting the number of timesteps from that\n",
"nt = 1000 # number of time steps after initial conditions\n", "nt = 3000 # number of time steps after initial conditions\n",
"\n", "\n",
"# derivatives of the pipeline constants\n", "# derivatives of the pipeline constants\n",
"dx = L/n # length of each pipe segment\n", "dx = L/n # length of each pipe segment\n",
"dt = dx/c # timestep according to method of characterisitics\n", "dt = dx/c # timestep according to method of characterisitics\n",
"nn = n+1 # number of nodes\n", "nn = n+1 # number of nodes\n",
"initial_level = 20. # water level in upstream reservoir [m]\n", "initial_level = 8. # water level in upstream reservoir [m]\n",
"p0 = rho*g*initial_level-v0**2*rho/2\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*dx,dx) # pl = pipe-length. position of the nodes on the pipeline\n",
"t_vec = np.arange(0,nt+1)*dt # time vector\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,n+1)*h_pipe/n # hydraulic head of pipeline at each node \n",
"v_init = np.full(nn,Q0/(A_pipe)) # initial velocity distribution in pipeline\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", "# 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",
"\n", "\n",
"\n", "\n",
"# reservoir\n", "# reservoir\n",
"# replace influx by time variable vector\n", "# replace influx by vector\n",
"initial_influx = Q0 # initial influx of volume to the reservoir [m³/s]\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_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_pipeline_pressure = p0 # Initial condition for the static pipeline pressure at the reservoir (= hydrostatic pressure - dynamic pressure) \n",
"initial_pressure_unit = 'Pa' # DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", "initial_pressure_unit = 'Pa' # DO NOT CHANGE! for pressure conversion in print statements and plot labels \n",
"conversion_pressure_unit = 'mWS' # 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 = 20. # total base are of the cuboid reservoir [m²] \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", "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_low = 0. # for yet-to-be-implemented warnings[m]\n",
"critical_level_high = np.inf # for yet-to-be-implemented warnings[m]\n", "critical_level_high = np.inf # for yet-to-be-implemented warnings[m]\n",
@@ -69,6 +74,7 @@
"# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\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 = 1000 # number of simulation steps of reservoir in between timesteps of pipeline \n",
"simulation_timestep = dt/nt_eRK4\n", "simulation_timestep = dt/nt_eRK4\n",
"\n",
"\n" "\n"
] ]
}, },
@@ -91,23 +97,24 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 13, "execution_count": 3,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# create objects\n", "# create objects\n",
"\n", "\n",
"V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,simulation_timestep)\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_steady_state(initial_influx,initial_level,initial_pressure_unit,conversion_pressure_unit)\n",
"V.set_influx(initial_influx)\n",
"V.set_outflux(initial_outflux)\n",
"V.set_pressure(initial_pipeline_pressure,initial_pressure_unit,conversion_pressure_unit)\n",
"\n", "\n",
"pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n", "pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n",
"pipe.set_pressure_propagation_velocity(c)\n", "pipe.set_pressure_propagation_velocity(c)\n",
"pipe.set_number_of_timesteps(nt)\n", "pipe.set_number_of_timesteps(nt)\n",
"pipe.set_initial_pressure(p_init,initial_pressure_unit,conversion_pressure_unit)\n", "pipe.set_steady_state(initial_influx,V.level,pl_vec,h_vec,initial_pressure_unit,conversion_pressure_unit)\n",
"pipe.set_initial_flow_velocity(v_init)\n", "\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", "\n",
"# display the attributes of the created reservoir and pipeline object\n", "# display the attributes of the created reservoir and pipeline object\n",
"# V.get_info(full=True)\n", "# V.get_info(full=True)\n",
@@ -116,15 +123,15 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 14, "execution_count": 4,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# initialization for timeloop\n", "# initialization for timeloop\n",
"\n", "\n",
"# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n", "# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n",
"v_old = v_init.copy()\n", "v_old = pipe.v0.copy()\n",
"p_old = p_init.copy()\n", "p_old = pipe.p0.copy()\n",
"\n", "\n",
"# prepare the vectors in which the temporal evolution of the boundary conditions are stored\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", " # keep in mind, that the velocity at the turbine and the pressure at the reservoir are set manually and\n",
@@ -134,36 +141,26 @@
"v_boundary_tur = 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_res = np.empty_like(t_vec)\n",
"p_boundary_tur = np.empty_like(t_vec)\n", "p_boundary_tur = np.empty_like(t_vec)\n",
"influx_vec = np.empty_like(t_vec)\n",
"\n", "\n",
"# prepare the vectors that store the temporal evolution of the level in the reservoir\n", "# prepare the vectors that store the temporal evolution of the level in the reservoir\n",
"level_vec = np.full(nt+1,initial_level) # level at the end of each pipeline timestep\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_2 = np.empty([nt_eRK4]) # level throughout each reservoir timestep-used for plotting and overwritten afterwards\n",
"\n", "\n",
"# set the boudary conditions for the first timestep\n", "# set the boundary conditions for the first timestep\n",
"v_boundary_res[0] = v_old[0]\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_res[0] = p_old[0]\n",
"p_boundary_tur[0] = p_old[-1]\n", "p_boundary_tur[0] = p_old[-1]\n",
"\n", "\n",
"v_boundary_tur[:] = v_old[0] \n", "LA_soll_vec = np.full_like(t_vec,T1.LA)\n",
"v_boundary_tur[1:] = 0 # instantaneous closing\n", "LA_soll_vec[1500:]= 0\n",
"\n", "\n",
"const = int(np.min([1000,round(nt/1.25)])) \n",
"# # v_boundary_tur[0:const] = np.linspace(v_old[-1],0,const) # linear closing\n",
"v_boundary_tur[0:const] = v_old[1]*np.cos(t_vec[0:const]*2*np.pi)**2 # oscillating\n",
"\n",
"influx_vec[0] = initial_influx # instantaneous closing\n",
"influx_vec[1:] = initial_influx # instantaneous closing\n",
"\n",
"const2 = int(np.min([1000,round(nt/1.25)])) \n",
"# influx_vec[0:const2] = np.linspace(v_old[-1],0,const2) # linear closing\n",
"# influx_vec[0:const2] = initial_influx*np.cos(t_vec[0:const2]*2*np.pi)**2 # oscillating\n",
"\n" "\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 15, "execution_count": 5,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -191,8 +188,15 @@
"# axs1[2].autoscale()\n", "# axs1[2].autoscale()\n",
"fig1.tight_layout()\n", "fig1.tight_layout()\n",
"fig1.show()\n", "fig1.show()\n",
"plt.pause(1)\n", "plt.pause(1)\n"
"\n", ]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# loop through time steps of the pipeline\n", "# loop through time steps of the pipeline\n",
"for it_pipe in range(1,pipe.nt+1):\n", "for it_pipe in range(1,pipe.nt+1):\n",
"\n", "\n",
@@ -200,7 +204,6 @@
" # set initial conditions for the reservoir time evolution calculted with e-RK4\n", " # set initial conditions for the reservoir time evolution calculted with e-RK4\n",
" V.pressure = p_old[0]\n", " V.pressure = p_old[0]\n",
" V.outflux_vel = v_old[0]\n", " V.outflux_vel = v_old[0]\n",
" V.influx = influx_vec[it_pipe]\n",
" # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\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", " for it_res in range(nt_eRK4):\n",
" V.e_RK_4() # call e-RK4 to update outflux\n", " V.e_RK_4() # call e-RK4 to update outflux\n",
@@ -208,14 +211,18 @@
" V.set_volume() # update volume in reservoir\n", " V.set_volume() # update volume in reservoir\n",
" level_vec_2[it_res] = V.level # save for plotting\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", " 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", " break \n",
" level_vec[it_pipe] = V.level \n", " level_vec[it_pipe] = V.level \n",
"\n", "\n",
" # set boundary conditions for the next timestep of the characteristic method\n", " # set boundary conditions for the next timestep of the characteristic method\n",
" p_boundary_res[it_pipe] = rho*g*V.level-v_old[1]**2*rho/2\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", " 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", " +dt*g*np.sin(alpha)\n",
"\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",
"\n",
" # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\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", " 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", " p_boundary_tur[it_pipe] = pipe.p_boundary_tur\n",
@@ -236,7 +243,7 @@
" fig1.canvas.draw()\n", " fig1.canvas.draw()\n",
" fig1.tight_layout()\n", " fig1.tight_layout()\n",
" fig1.show()\n", " fig1.show()\n",
" plt.pause(0.00001) \n", " plt.pause(0.1) \n",
"\n", "\n",
" # prepare for next loop\n", " # prepare for next loop\n",
" p_old = pipe.p_old\n", " p_old = pipe.p_old\n",
@@ -248,7 +255,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 16, "execution_count": 7,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [

324
Untertweng_mit_Pegelregler.ipynb Normal file
View File

@@ -0,0 +1,324 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from functions.pressure_conversion import pressure_conversion\n",
"from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class\n",
"from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class\n",
"from Turbinen.Turbinen_class_file import Francis_Turbine\n",
"from Regler.Regler_class_file import PI_controller_class"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"#define constants\n",
"\n",
"#Turbine\n",
"Q_nenn = 0.85\n",
"p_nenn,_ = pressure_conversion(10.6,'bar','Pa')\n",
"\n",
"# physics\n",
"g = 9.81 # gravitational acceleration [m/s²]\n",
"rho = 1000. # density of water [kg/m³]\n",
"\n",
"# pipeline\n",
"L = 535.+478. # length of pipeline [m]\n",
"D = 0.9 # pipe diameter [m]\n",
"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",
"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",
"\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",
"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",
"\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_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",
"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",
"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,
"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",
"\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",
"\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",
"\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()"
]
},
{
"cell_type": "code",
"execution_count": 25,
"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",
"\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.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",
"\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",
"\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",
"\n",
"LA_soll_vec = np.full_like(t_vec,T1.LA)\n",
"Pegelregler.control_variable = T1.LA\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib qt5\n",
"# time loop\n",
"\n",
"# create a figure and subplots to display the velocity and pressure distribution across the pipeline in each pipeline step\n",
"fig1,axs1 = plt.subplots(2,1)\n",
"fig1.suptitle(str(0) +' s / '+str(round(t_vec[-1],2)) + ' s' )\n",
"axs1[0].set_title('Pressure distribution in pipeline')\n",
"axs1[1].set_title('Velocity distribution in pipeline')\n",
"axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"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",
"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",
"fig1.tight_layout()\n",
"fig1.show()\n",
"plt.pause(1)\n"
]
},
{
"cell_type": "code",
"execution_count": 27,
"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",
"\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",
" # 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",
"\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",
"\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",
"\n",
" # perform the next timestep via the characteristic method\n",
" pipe.timestep_characteristic_method()\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_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",
"\n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": 28,
"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].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[0,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
"axs2[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
"axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[1,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
"axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\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",
"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
}

22
untertweng.txt Normal file
View File

@@ -0,0 +1,22 @@
L = 535 m dn 800 mm
478 m dn 1000 mm
Ersatzdurchmesser
h_pipe
h 851.78 Pegel + Leitungsgefälle
Leitungsgefälle: 113
Fläche 4.25x10.5 + 30m² = 74 m²
Pegelminimum: 851.18 m
Unterwasserpegel 738.56
Gesamtfallhöhe = 851.78-738.56
Rohrreibung: 0.014 f_D = lambda
c = 500 m/s
Q_0 = 100%*0.75+30%*0.75
Q_extrem = 30%*0.75
Q = LA*Q_nenn*sqrt(H/H_n)