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

code cleanup:

Browse Source 
consistenly use getter and setter methods
commenting etc
This commit is contained in:
Brantegger Georg
2022-07-27 11:40:58 +02:00
parent ac8bfdb7c6
commit d1c15090dc
13 changed files with 956 additions and 584 deletions
Download Patch File Download Diff File Expand all files Collapse all files

33
Ausgleichsbecken/Ausgleichsbecken_class_file.py
View File

@@ -44,8 +44,8 @@ class Ausgleichsbecken_class:
density_unit_print = 'kg/m³'
flux_unit_print = 'm³/s'
level_unit_print = 'm'
time_unit_print = 's'
pressure_unit_print = '--' # will be set by .set_pressure() method
time_unit_print = 's'
velocity_unit_print = 'm/s'
volume_unit_print = 'm³'
@@ -102,8 +102,7 @@ class Ausgleichsbecken_class:
# 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_outflux_vel = ss_outflux/self.area_outflux
ss_pressure = self.density*self.g*ss_level-ss_outflux_vel**2*self.density/2
ss_pressure = self.density*self.g*ss_level-(ss_outflux/self.area_outflux)**2*self.density/2
self.set_influx(ss_influx)
self.set_initial_level(ss_level)
@@ -169,6 +168,9 @@ class Ausgleichsbecken_class:
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
@@ -178,31 +180,32 @@ class Ausgleichsbecken_class:
# sets volume in reservoir based on self.level
return self.level*self.area
def update_pressure(self):
p_new = self.density*self.g*self.level-(self.outflux/self.area_outflux)**2*self.density/2
return p_new
def timestep_reservoir_evolution(self):
# update outflux and outflux velocity based on current pipeline pressure and waterlevel in reservoir
yn = self.outflux/self.area_outflux # outflux velocity
h = self.level
dt = self.timestep
p = self.pressure
# assume constant pipeline pressure during timestep
# e_RK_4 timestep is way smalle than timestep of characteristic method, so this should be a valid approx.
# (furthermore I have no idea how to approximate p_hs otherwise :/ )
p_hs = self.pressure
A_a = self.area_outflux
A = self.area
h_hs = self.update_level(dt/2)
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,A,A_a,self.pressure,rho,g)
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 = ynp1*self.area_outflux
self.outflux = ynp1*A_a
self.level = self.update_level(dt)
self.volume = self.update_volume()
self.pressure = self.update_pressure()

36
Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb
View File

@@ -28,9 +28,9 @@
"# 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",
"# 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",
@@ -46,7 +46,7 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
@@ -72,7 +72,6 @@
"i_max = -1\n",
"\n",
"for i in range(np.size(time_vec)-1):\n",
" # update to include p_halfstep\n",
" V.set_pressure(pressure_vec[i])\n",
" V.timestep_reservoir_evolution()\n",
" outflux_vec[i+1] = V.get_current_outflux()\n",
@@ -85,12 +84,12 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"\n",
"fig1, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1)\n",
"fig1, (ax1, ax2, ax3) = plt.subplots(3, 1)\n",
"fig1.set_figheight(10)\n",
"fig1.suptitle('Ausgleichsbecken')\n",
"\n",
@@ -109,29 +108,6 @@
"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() "
]

21
Druckrohrleitung/Druckrohrleitung_class_file.py
View File

@@ -33,6 +33,8 @@ class Druckrohrleitung_class:
self.density = rho # density of the liquid in the pipeline
self.g = g # gravitational acceleration
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
@@ -98,23 +100,22 @@ class Druckrohrleitung_class:
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
self.v_boundary_tur = v_turbine # at new timestep
self.p_boundary_res = p_reservoir # at new timestep
v_boundary_tur = v_turbine # at new timestep
p_boundary_res = p_reservoir # at new timestep
# calculate the missing boundary conditions
self.v_boundary_res = v_old_res+1/(rho*c)*(p_reservoir-p_old_res)+dt*g*np.sin(alpha)-f_D*dt/(2*D)*abs(v_old_res)*v_old_res
self.p_boundary_tur = p_old_tur-rho*c*(v_turbine-v_old_tur)+rho*c*dt*g*np.sin(alpha)-f_D*rho*c*dt/(2*D)*abs(v_old_tur)*v_old_tur
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
# write boundary conditions to the velocity/pressure vectors of the next timestep
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()
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,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.dia**2/4*np.pi))
ss_v0 = np.full(self.n_seg+1,ss_flux/self.A)
# the static pressure is given by the hydrostatic pressure, corrected for friction losses and dynamic pressure
ss_pressure = (self.density*self.g*(ss_level_reservoir+h_vec)-ss_v0**2*self.density/2)-(self.f_D*pl_vec/self.dia*self.density/2*ss_v0**2)

230
Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
@@ -16,64 +16,98 @@
"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"
"from functions.pressure_conversion import pressure_conversion\n",
"from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class"
]
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib qt5\n",
"#define constants\n",
"#define constants pipe\n",
"\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",
"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",
"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 = 5000 # 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",
"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_influx = 2. # m³/s\n",
"initial_level = 10. # m\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*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"
"h_vec = np.arange(0,h_pipe+h_pipe/n,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,
"execution_count": null,
"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_influx,initial_level,pl_vec,h_vec)"
"pipe.set_steady_state(initial_flux,initial_level,pl_vec,h_vec)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(V.get_current_influx())\n",
"print(V.get_current_outflux())\n",
"print(V.get_current_level())\n",
"print(V.get_current_pressure())\n",
"print(pipe.get_current_pressure_distribution()[0])\n",
"print(pipe.get_current_velocity_distribution()*area_pipe)\n",
"print(pipe.get_current_velocity_distribution())"
]
},
{
"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",
@@ -94,80 +128,110 @@
"p_boundary_tur[0] = p_old[-1]\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"fig2,axs2 = plt.subplots(2,1)\n",
"axs2[0].set_title('Pressure distribution in pipeline')\n",
"axs2[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs2[0].set_ylabel(r'$p$ [mWS]')\n",
"lo_00, = axs2[0].plot(pl_vec,pressure_conversion(p_old,'Pa','mWS'),marker='.')\n",
"axs2[0].set_ylim([0.9*np.min(pressure_conversion(p_old,'Pa','mWS')),1.1*np.max(pressure_conversion(p_old,'Pa','mWS'))])\n",
"\n",
"axs2[1].set_title('Velocity distribution in pipeline')\n",
"axs2[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs2[1].set_ylabel(r'$p$ [mWS]')\n",
"lo_01, = axs2[1].plot(pl_vec,v_old,marker='.')\n",
"axs2[1].set_ylim([0.9*np.min(v_old),1.1*np.max(v_boundary_res)])\n",
"\n",
"fig2.tight_layout()\n",
"plt.pause(5)\n",
"\n",
"\n",
"for it in range(1,pipe.nt):\n",
" pipe.set_boundary_conditions_next_timestep(p_boundary_res[0],v_boundary_tur[0])\n",
" pipe.timestep_characteristic_method()\n",
" lo_00.set_ydata(pressure_conversion(pipe.get_current_pressure_distribution(),'Pa','mWS'))\n",
" lo_01.set_ydata(pipe.get_current_velocity_distribution())\n",
"\n",
" v_boundary_res[it] = pipe.get_current_velocity_distribution()[0]\n",
" v_boundary_tur[it] = pipe.get_current_velocity_distribution()[-1]\n",
" p_boundary_res[it] = pipe.get_current_pressure_distribution()[0]\n",
" p_boundary_tur[it] = pipe.get_current_pressure_distribution()[-1]\n",
"\n",
"\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": 6,
"metadata": {},
"outputs": [],
"source": [
"fig3,axs3 = plt.subplots(2,2)\n",
"axs3[0,0].set_title('Pressure Reservoir')\n",
"axs3[0,0].plot(t_vec,pressure_conversion(p_boundary_res,'Pa','mWS'))\n",
"axs3[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs3[0,0].set_ylabel(r'$p$ [mWS]')\n",
"axs3[0,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_res,'Pa','mWS')),1.1*np.max(pressure_conversion(p_boundary_res,'Pa','mWS'))])\n",
"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",
"axs3[0,1].set_title('Velocity Reservoir')\n",
"axs3[0,1].plot(t_vec,v_boundary_res)\n",
"axs3[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs3[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
"axs3[0,1].set_ylim([0.9*np.min(v_boundary_res),1.1*np.max(v_boundary_res)])\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",
"axs3[1,0].set_title('Pressure Turbine')\n",
"axs3[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,'Pa','mWS'))\n",
"axs3[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs3[1,0].set_ylabel(r'$p$ [mWS]')\n",
"axs3[1,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_tur,'Pa','mWS')),1.1*np.max(pressure_conversion(p_boundary_tur,'Pa','mWS'))])\n",
"fig1.tight_layout()\n",
"plt.pause(1)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"\n",
"axs3[1,1].set_title('Velocity Turbine')\n",
"axs3[1,1].plot(t_vec,v_boundary_tur)\n",
"axs3[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs3[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
"axs3[1,1].set_ylim([0.9*np.min(v_boundary_tur),1.1*np.max(v_boundary_tur)])\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",
"fig3.tight_layout()\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(t_vec[-1]))\n",
" fig1.canvas.draw()\n",
" fig1.tight_layout()\n",
" plt.pause(0.00001)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"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()"
]
}

19
Kraftwerk/Kraftwerk_class_file.py Normal file
View 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
View 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
}

275
Regler/regler_test.ipynb → Regler/Pegelregler_test.ipynb
View File

@@ -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.e_RK_4() \n",
" V.level = V.update_level(V.timestep) \n",
" V.set_volume() \n",
" level_vec[i] = V.level \n",
" V.set_outflux(Q_vec[i])\n",
"\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())"
]
}
],

111
Regler/Regler_class_file.py
View File

@@ -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,58 +74,47 @@ 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
# 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
# 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:
@@ -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)

121
Turbinen/Turbinen_class_file.py
View File

@@ -1,3 +1,4 @@
from time import time
import numpy as np
#importing pressure conversion function
import sys
@@ -8,35 +9,117 @@ 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 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
def get_Q(self):
#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
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
self.set_LA(self.LA-LA_diff,display_warning=False) # set new LA
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)

78
Turbinen/Turbinen_test.ipynb → Turbinen/Turbinen_test_Kennlinie.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
@@ -23,18 +23,39 @@
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 5,
"metadata": {},
"outputs": [],
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"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 = -1 s \n",
"Current flux = -1 m³/s \n",
"Current pipe pressure = -1.0 mWS \n",
"Current LA = -1.0 % \n",
"Simulation timestep = -1 s \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)"
"Untertweng1 = Francis_Turbine(Q_nenn,p_nenn)\n",
"Untertweng1.get_info(full=True)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
@@ -47,7 +68,7 @@
},
{
"cell_type": "code",
"execution_count": 20,
"execution_count": 7,
"metadata": {},
"outputs": [
{
@@ -56,14 +77,14 @@
"Text(0.5, 0, 'Q [m³/s]')"
]
},
"execution_count": 20,
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "7650cd2a1f9a4b87b98e70add29f11a7",
"model_id": "75adb3cb47e642e3a5606cb41efedf72",
"version_major": 2,
"version_minor": 0
},
@@ -94,8 +115,8 @@
"for i in range(n_p):\n",
" for j in range(n_LA):\n",
" Untertweng1.set_pressure(pp[i,j])\n",
" Untertweng1.set_LA(ll[i,j])\n",
" Q_mat[i,j] = Untertweng1.get_Q()\n",
" Untertweng1.set_LA(ll[i,j],display_warning=False)\n",
" Q_mat[i,j] = Untertweng1.get_current_Q()\n",
"\n",
"fig1 = plt.figure()\n",
"ax1 = plt.axes(projection='3d')\n",
@@ -108,31 +129,23 @@
},
{
"cell_type": "code",
"execution_count": 27,
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\BRANT\\AppData\\Local\\Temp\\9\\ipykernel_7508\\1599598770.py:5: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
" fig = plt.figure()\n"
]
},
{
"data": {
"text/plain": [
"(0.0, 1.275)"
]
},
"execution_count": 27,
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "52152d0b96c74a4ebfb18041a22d8d0e",
"model_id": "224f00f9bf85446b845685a08ed27c68",
"version_major": 2,
"version_minor": 0
},
@@ -170,7 +183,7 @@
},
{
"cell_type": "code",
"execution_count": 30,
"execution_count": 9,
"metadata": {},
"outputs": [
{
@@ -179,14 +192,14 @@
"(0.0, 1.275)"
]
},
"execution_count": 30,
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "bf6b0fe84d264693813f6e991600ece1",
"model_id": "92741b823c9749c9820ee7b5ba47a6bc",
"version_major": 2,
"version_minor": 0
},
@@ -221,23 +234,6 @@
"plt.title('P = '+ str(p_test2) + ' [Pa]')\n",
"plt.ylim([0,1.5*Q_nenn])"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"53000.0\n"
]
}
],
"source": [
"print(pp[10,5])"
]
}
],
"metadata": {

90
Turbinen/messy.ipynb Normal file
View 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
}

166
Untertweng.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
@@ -17,15 +17,16 @@
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#define constants\n",
"\n",
"#Turbine\n",
"Q_nenn = 0.85\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",
@@ -37,33 +38,29 @@
"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",
"n = 200 # 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",
"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 = 20000 # 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",
"area_base = 74. # total base are of the cuboid reservoir [m²] \n",
@@ -78,43 +75,28 @@
"\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,
"execution_count": null,
"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,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)\n",
"pipe.set_steady_state(initial_flux,initial_level,pl_vec,h_vec)\n",
"\n",
"initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n",
"\n",
"T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time,timestep = dt)\n",
"T1.set_steady_state(initial_flux,initial_pressure_turbine)\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",
@@ -123,19 +105,19 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": null,
"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.zeros_like(t_vec)\n",
"v_boundary_tur = np.zeros_like(t_vec)\n",
@@ -143,8 +125,7 @@
"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.zeros([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,14 +134,15 @@
"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"
]
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
@@ -176,16 +158,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),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"
@@ -193,60 +170,62 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# loop through time steps of the pipeline\n",
"for it_pipe in range(1,pipe.nt+1):\n",
"\n",
" if t_vec[it_pipe]>20:\n",
" if V.get_current_influx() > 0:\n",
" V.set_influx(np.max([V.get_current_influx()-initial_flux*5*1e-3,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.update_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",
"\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(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n",
" p_boundary_tur[it_pipe] = pipe.p_boundary_tur\n",
" v_boundary_res[it_pipe] = pipe.v_boundary_res\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),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",
" "
@@ -254,35 +233,44 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": null,
"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))\n",
"axs2[0,1].plot(t_vec,v_boundary_res)\n",
"axs2[0,1].set_ylim(-2*Q_nenn,+2*Q_nenn)\n",
"axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,initial_pressure_unit, conversion_pressure_unit))\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()"
]

235
Untertweng_mit_Pegelregler.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
"execution_count": 22,
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
@@ -18,20 +18,29 @@
},
{
"cell_type": "code",
"execution_count": 23,
"execution_count": 19,
"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 = 70 #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 = 10.\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,20 +51,22 @@
"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 = 1000 # 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",
"area_base = 74. # total base are of the cuboid reservoir [m²] \n",
@@ -63,13 +74,6 @@
"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",
@@ -77,75 +81,55 @@
"\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": 20,
"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,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",
"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": 22,
"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",
"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,15 +138,14 @@
"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",
"LA_ist_vec = np.full_like(t_vec,T1.LA)\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 26,
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
@@ -178,16 +161,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"
@@ -195,65 +173,83 @@
},
{
"cell_type": "code",
"execution_count": 27,
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.2146505196687856\n"
]
}
],
"source": [
"print(initial_flux/area_outflux)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"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",
" if it_pipe > 0.015*(nt+1):\n",
" if V.get_current_influx() > 0:\n",
" V.set_influx(np.max([V.get_current_influx()-initial_flux*5*1e-3,0.]))\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]\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",
" # 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",
" 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",
" 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",
" # 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,37 +257,64 @@
},
{
"cell_type": "code",
"execution_count": 28,
"execution_count": 26,
"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": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0\n"
]
}
],
"source": [
"print(np.sum(v_boundary_res[2500:])*area_outflux)"
]
}
],
"metadata": {