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

further code cleanup

Browse Source
This commit is contained in:
Brantegger Georg
2022-07-25 15:59:46 +02:00
parent 0de946f8ac
commit ac8bfdb7c6
9 changed files with 448 additions and 438 deletions
Download Patch File Download Diff File Expand all files Collapse all files

37
Ausgleichsbecken/Ausgleichsbecken_class_file.py
View File

@@ -69,15 +69,11 @@ class Ausgleichsbecken_class:
# setter # setter
def update_volume(self):
# sets volume in reservoir based on self.level
self.volume = self.level*self.area
def set_initial_level(self,initial_level): def set_initial_level(self,initial_level):
# sets the level in the reservoir and should only be called during initialization # sets the level in the reservoir and should only be called during initialization
if self.level == '--': if self.level == '--':
self.level = initial_level self.level = initial_level
self.update_volume() self.volume = self.update_volume()
else: else:
raise Exception('Initial level was already set once. Use the .update_level(self,timestep) method to update level based on net flux.') raise Exception('Initial level was already set once. Use the .update_level(self,timestep) method to update level based on net flux.')
@@ -91,28 +87,33 @@ class Ausgleichsbecken_class:
# positive outflux means that liquid flows out of reservoir the reservoir # positive outflux means that liquid flows out of reservoir the reservoir
self.outflux = outflux self.outflux = outflux
def set_pressure(self,pressure,display_pressure_unit): def set_initial_pressure(self,pressure,display_pressure_unit):
# sets the static pressure present at the outlet of the reservoir # sets the static pressure present at the outlet of the reservoir
# units are used to convert and display the pressure # units are used to convert and display the pressure
self.pressure = pressure self.pressure = pressure
self.pressure_unit_print = display_pressure_unit self.pressure_unit_print = display_pressure_unit
def set_pressure(self,pressure):
# sets the static pressure present at the outlet of the reservoir
# units are used to convert and display the pressure
self.pressure = pressure
def set_steady_state(self,ss_influx,ss_level,display_pressure_unit): def set_steady_state(self,ss_influx,ss_level,display_pressure_unit):
# find the steady state (ss) condition in which the net flux is zero # set the steady state (ss) condition in which the net flux is zero
# set pressure acting on the outflux so that the level stays constant # set pressure acting on the outflux area so that the level stays constant
ss_outflux = ss_influx ss_outflux = ss_influx
ss_outflux_vel = ss_outflux/self.area_outflux 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_vel**2*self.density/2
self.set_initial_level(ss_level)
self.set_influx(ss_influx) self.set_influx(ss_influx)
self.set_initial_level(ss_level)
self.set_initial_pressure(ss_pressure,display_pressure_unit)
self.set_outflux(ss_outflux) self.set_outflux(ss_outflux)
self.set_pressure(ss_pressure,display_pressure_unit)
# getter # getter
def get_info(self, full = False): def get_info(self, full = False):
new_line = '\n' new_line = '\n'
p,_ = pressure_conversion(self.pressure,self.pressure_unit,self.pressure_unit_print) p = pressure_conversion(self.pressure,self.pressure_unit,self.pressure_unit_print)
if full == True: if full == True:
@@ -164,6 +165,7 @@ class Ausgleichsbecken_class:
# methods # methods
def update_level(self,timestep): def update_level(self,timestep):
# update level based on net flux and timestep by calculating the volume change in # 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 # the timestep and the converting the new volume to a level by assuming a cuboid reservoir
@@ -172,8 +174,12 @@ class Ausgleichsbecken_class:
new_level = (self.volume+delta_V)/self.area new_level = (self.volume+delta_V)/self.area
return new_level return new_level
def update_volume(self):
# sets volume in reservoir based on self.level
return self.level*self.area
def e_RK_4(self):
def timestep_reservoir_evolution(self):
# update outflux and outflux velocity based on current pipeline pressure and waterlevel in reservoir # update outflux and outflux velocity based on current pipeline pressure and waterlevel in reservoir
yn = self.outflux/self.area_outflux # outflux velocity yn = self.outflux/self.area_outflux # outflux velocity
h = self.level h = self.level
@@ -196,4 +202,7 @@ class Ausgleichsbecken_class:
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)+ \ 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)) 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*self.area_outflux
self.level = self.update_level(dt)
self.volume = self.update_volume()

30
Ausgleichsbecken/Ausgleichsbecken_test.ipynb → Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 4, "execution_count": 1,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -21,7 +21,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 5, "execution_count": 2,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -46,7 +46,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 6, "execution_count": 4,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -56,29 +56,27 @@
"# 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",
"# converted_pressure,_ = pressure_conversion(initial_pipeline_pressure,input_unit = initial_pressure_unit, target_unit = 'Pa')\n", "# V.set_initial_pressure(pressure_conversion(initial_pipeline_pressure,input_unit = initial_pressure_unit, target_unit = 'Pa'),conversion_pressure_unit)\n",
"# V.pressure = converted_pressure\n", "# V.pressure = converted_pressure\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", "\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] = V.outflux\n", "outflux_vec[0] = V.get_current_outflux()\n",
"level_vec = np.empty_like(time_vec)\n", "level_vec = np.empty_like(time_vec)\n",
"level_vec[0] = V.level\n", "level_vec[0] = V.get_current_level()\n",
"\n", "\n",
"# pressure_vec = np.full_like(time_vec,converted_pressure)*((np.sin(time_vec)+1)*np.exp(-time_vec/50))\n", "# pressure_vec = np.full_like(time_vec,converted_pressure)*((np.sin(time_vec)+1)*np.exp(-time_vec/50))\n",
"pressure_vec = np.full_like(time_vec,V.pressure)\n", "pressure_vec = np.full_like(time_vec,V.get_current_pressure())\n",
" \n", " \n",
"i_max = -1\n", "i_max = -1\n",
"\n", "\n",
"for i in range(np.size(time_vec)-1):\n", "for i in range(np.size(time_vec)-1):\n",
" # update to include p_halfstep\n", " # update to include p_halfstep\n",
" V.pressure = pressure_vec[i]\n", " V.set_pressure(pressure_vec[i])\n",
" V.e_RK_4()\n", " V.timestep_reservoir_evolution()\n",
" V.level = V.update_level(V.timestep)\n", " outflux_vec[i+1] = V.get_current_outflux()\n",
" V.update_volume()\n", " level_vec[i+1] = V.get_current_level()\n",
" outflux_vec[i+1] = V.outflux\n",
" level_vec[i+1] = V.level\n",
" 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",
@@ -87,7 +85,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 7, "execution_count": 8,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -106,7 +104,7 @@
"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_conversion(pressure_vec[:i_max],'Pa',conversion_pressure_unit)[0], label='Pipeline pressure at reservoir')\n", "ax3.plot(time_vec[:i_max],pressure_conversion(pressure_vec[:i_max],'Pa',conversion_pressure_unit), label='Pipeline pressure at reservoir')\n",
"ax3.set_ylabel(r'$p_{pipeline}$ ['+conversion_pressure_unit+']')\n", "ax3.set_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",

165
Druckrohrleitung/Druckrohrleitung_class_file.py
View File

@@ -1,51 +1,40 @@
import numpy as np import numpy as np
#importing pressure conversion function
import sys
import os
current = os.path.dirname(os.path.realpath(__file__))
parent = os.path.dirname(current)
sys.path.append(parent)
from functions.pressure_conversion import pressure_conversion
class Druckrohrleitung_class: class Druckrohrleitung_class:
# units # units
acceleration_unit = r'$\mathrm{m}/\mathrm{s}^2$' acceleration_unit = r'$\mathrm{m}/\mathrm{s}^2$'
angle_unit = 'rad' angle_unit = 'rad'
area_unit = r'$\mathrm{m}^2$' area_unit = r'$\mathrm{m}^2$'
density_unit = r'$\mathrm{kg}/\mathrm{m}^3$' density_unit = r'$\mathrm{kg}/\mathrm{m}^3$'
flux_unit = r'$\mathrm{m}^3/\mathrm{s}$' flux_unit = r'$\mathrm{m}^3/\mathrm{s}$'
length_unit = 'm' length_unit = 'm'
time_unit = 's' pressure_unit = 'Pa'
velocity_unit = r'$\mathrm{m}/\mathrm{s}$' # for flux and pressure propagation time_unit = 's'
volume_unit = r'$\mathrm{m}^3$' velocity_unit = r'$\mathrm{m}/\mathrm{s}$' # for flux and pressure propagation
volume_unit = r'$\mathrm{m}^3$'
acceleration_unit_print = 'm/s²' acceleration_unit_print = 'm/s²'
angle_unit_print = 'rad' angle_unit_print = 'rad'
area_unit_print = 'm²' area_unit_print = 'm²'
density_unit_print = 'kg/m³' density_unit_print = 'kg/m³'
flux_unit_print = 'm³/s' flux_unit_print = 'm³/s'
length_unit_print = 'm' length_unit_print = 'm'
time_unit_print = 's' time_unit_print = 's'
velocity_unit_print = 'm/s' # for flux and pressure propagation velocity_unit_print = 'm/s' # for flux and pressure propagation
volume_unit_print = 'm³' volume_unit_print = 'm³'
# init # init
def __init__(self,total_length,diameter,number_segments,pipeline_angle,Darcy_friction_factor,rho=1000,g=9.81): def __init__(self,total_length,diameter,number_segments,pipeline_angle,Darcy_friction_factor,rho=1000,g=9.81):
self.length = total_length self.length = total_length # total length of the pipeline
self.dia = diameter self.dia = diameter # diameter of the pipeline
self.n_seg = number_segments self.n_seg = number_segments # number of segments for the method of characteristics
self.angle = pipeline_angle self.angle = pipeline_angle # angle of the pipeline
self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient
self.rho = rho self.density = rho # density of the liquid in the pipeline
self.g = g self.g = g # gravitational acceleration
self.dx = total_length/number_segments self.dx = total_length/number_segments # length of each segment
self.l_vec = np.arange(0,(number_segments+1)*self.dx,self.dx) self.l_vec = np.arange(0,(number_segments+1),1)*self.dx # vector giving the distance from each node to the start of the pipeline
# initialize for get_info method # initialize for get_info method
self.c = '--' self.c = '--'
@@ -53,32 +42,32 @@ class Druckrohrleitung_class:
# setter # setter
def set_pressure_propagation_velocity(self,c): def set_pressure_propagation_velocity(self,c):
self.c = c self.c = c # propagation velocity of the pressure wave
self.dt = self.dx/c self.dt = self.dx/c # timestep derived from c, demanded by the method of characteristics
def set_number_of_timesteps(self,number_timesteps): def set_number_of_timesteps(self,number_timesteps):
self.nt = number_timesteps self.nt = number_timesteps # number of timesteps
if self.c == '--': if self.c == '--':
raise Exception('Please set the pressure propagation velocity before setting the number of timesteps.') raise Exception('Please set the pressure propagation velocity before setting the number of timesteps.')
else: else:
self.t_vec = np.arange(0,self.nt*self.dt,self.dt) self.t_vec = np.arange(0,self.nt*self.dt,self.dt)
def set_initial_pressure(self,pressure,pressure_unit,display_pressure_unit): def set_initial_pressure(self,pressure):
if np.size(pressure) == 1: # initialize the pressure distribution in the pipeline
if np.size(pressure) == 1:
self.p0 = np.full_like(self.l_vec,pressure) self.p0 = np.full_like(self.l_vec,pressure)
elif np.size(pressure) == np.size(self.l_vec): elif np.size(pressure) == np.size(self.l_vec):
self.p0 = pressure self.p0 = pressure
else: else:
raise Exception('Unable to assign initial pressure. Input has to be of size 1 or' + np.size(self.l_vec)) raise Exception('Unable to assign initial pressure. Input has to be of size 1 or' + np.size(self.l_vec))
self.pressure_unit = pressure_unit
self.pressure_unit_print = display_pressure_unit
#initialize the vectors in which the old and new pressures are stored for the method of characteristics #initialize the vectors in which the old and new pressures are stored for the method of characteristics
self.p_old = self.p0.copy() self.p_old = self.p0.copy()
self.p = np.empty_like(self.p_old) self.p = self.p0.copy()
def set_initial_flow_velocity(self,velocity): def set_initial_flow_velocity(self,velocity):
if np.size(velocity) == 1: # initialize the velocity distribution in the pipeline
if np.size(velocity) == 1:
self.v0 = np.full_like(self.l_vec,velocity) self.v0 = np.full_like(self.l_vec,velocity)
elif np.size(velocity) == np.size(self.l_vec): elif np.size(velocity) == np.size(self.l_vec):
self.v0 = velocity self.v0 = velocity
@@ -86,35 +75,51 @@ class Druckrohrleitung_class:
raise Exception('Unable to assign initial velocity. Input has to be of size 1 or' + np.size(self.l_vec)) raise Exception('Unable to assign initial velocity. Input has to be of size 1 or' + np.size(self.l_vec))
#initialize the vectors in which the old and new velocities are stored for the method of characteristics #initialize the vectors in which the old and new velocities are stored for the method of characteristics
self.v_old = self.v0.copy() self.v_old = self.v0.copy()
self.v = np.empty_like(self.v_old) self.v = self.v0.copy()
def set_boundary_conditions_next_timestep(self,v_reservoir,p_reservoir,v_turbine): def set_boundary_conditions_next_timestep(self,p_reservoir,v_turbine):
rho = self.rho # derived from the method of characteristics, one can set the boundary conditions for the pressures and flow velocities at the reservoir and the turbine
# the boundary velocity at the turbine is specified by the flux through the turbine or an external boundary condition
# the pressure at the turbine will be calculated using the forward characteristic
# the boundary pressure at the reservoir is specified by the level in the reservoir of an external boundary condition
# the velocity at the reservoir will be calculated using the backward characteristic
# constants for a cleaner formula
rho = self.density
c = self.c c = self.c
f_D = self.f_D f_D = self.f_D
dt = self.dt dt = self.dt
D = self.dia D = self.dia
g = self.g g = self.g
alpha = self.angle alpha = self.angle
p_old = self.p_old[-2] # @ second to last node (the one before the turbine) p_old_tur = self.p_old[-2] # @ second to last node (the one before the turbine)
v_old = self.v_old[-2] # @ second to last node (the one before the turbine) v_old_tur = self.v_old[-2] # @ second to last node (the one before the turbine)
self.v_boundary_res = v_reservoir # at new timestep 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.v_boundary_tur = v_turbine # at new timestep
self.p_boundary_res = p_reservoir self.p_boundary_res = p_reservoir # at new timestep
self.p_boundary_tur = p_old-rho*c*(v_turbine-v_old)+rho*c*dt*g*np.sin(alpha)-f_D*rho*c*dt/(2*D)*abs(v_old)*v_old # 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
# write boundary conditions to the velocity/pressure vectors of the next timestep
self.v[0] = self.v_boundary_res.copy() self.v[0] = self.v_boundary_res.copy()
self.v[-1] = self.v_boundary_tur.copy() self.v[-1] = self.v_boundary_tur.copy()
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): 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.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) # 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)
self.set_initial_flow_velocity(ss_v0) self.set_initial_flow_velocity(ss_v0)
self.set_initial_pressure(ss_pressure,pressure_unit,display_pressure_unit) self.set_initial_pressure(ss_pressure)
# getter # getter
def get_info(self): def get_info(self):
@@ -142,30 +147,27 @@ class Druckrohrleitung_class:
f"Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object") f"Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object")
print(print_str) print(print_str)
def get_boundary_conditions_next_timestep(self): def get_current_pressure_distribution(self):
print('The pressure at the reservoir for the next timestep is', '\n', \ return self.p
pressure_conversion(self.p_boundary_res,self.pressure_unit,self.pressure_unit_print), '\n', \
'The velocity at the reservoir for the next timestep is', '\n', \
self.v_boundary_res, self.velocity_unit_print, '\n', \
'The pressure at the turbine for the next timestep is', '\n', \
pressure_conversion(self.p_boundary_tur,self.pressure_unit,self.pressure_unit_print), '\n', \
'The velocity at the turbine for the next timestep is', '\n', \
self.v_boundary_tur, self.velocity_unit_print)
def get_current_velocity_distribution(self):
return self.v
def timestep_characteristic_method(self): def timestep_characteristic_method(self):
#number of nodes # use the method of characteristics to calculate the pressure and velocities at all nodes except the boundary ones
nn = self.n_seg+1 # they are set with the .set_boundary_conditions_next_timestep() method beforehand
rho = self.rho
c = self.c nn = self.n_seg+1 # number of nodes
f_D = self.f_D rho = self.density # density of liquid
dt = self.dt c = self.c # pressure propagation velocity
D = self.dia f_D = self.f_D # Darcy friction coefficient
g = self.g dt = self.dt # timestep
alpha = self.angle D = self.dia # pipeline diametet
g = self.g # graviational acceleration
alpha = self.angle # pipeline angle
# Vectorize this loop?
for i in range(1,nn-1): for i in range(1,nn-1):
self.v[i] = 0.5*(self.v_old[i+1]+self.v_old[i-1])-0.5/(rho*c)*(self.p_old[i+1]-self.p_old[i-1]) \ self.v[i] = 0.5*(self.v_old[i+1]+self.v_old[i-1])-0.5/(rho*c)*(self.p_old[i+1]-self.p_old[i-1]) \
+dt*g*np.sin(alpha)-f_D*dt/(4*D)*(abs(self.v_old[i+1])*self.v_old[i+1]+abs(self.v_old[i-1])*self.v_old[i-1]) +dt*g*np.sin(alpha)-f_D*dt/(4*D)*(abs(self.v_old[i+1])*self.v_old[i+1]+abs(self.v_old[i-1])*self.v_old[i-1])
@@ -173,5 +175,8 @@ class Druckrohrleitung_class:
self.p[i] = 0.5*(self.p_old[i+1]+self.p_old[i-1]) - 0.5*rho*c*(self.v_old[i+1]-self.v_old[i-1]) \ self.p[i] = 0.5*(self.p_old[i+1]+self.p_old[i-1]) - 0.5*rho*c*(self.v_old[i+1]-self.v_old[i-1]) \
+f_D*rho*c*dt/(4*D)*(abs(self.v_old[i+1])*self.v_old[i+1]-abs(self.v_old[i-1])*self.v_old[i-1]) +f_D*rho*c*dt/(4*D)*(abs(self.v_old[i+1])*self.v_old[i+1]-abs(self.v_old[i-1])*self.v_old[i-1])
# prepare for next call
# use .copy() to write data to another memory location and avoid the usual python reference pointer
# else one can overwrite data by accidient and change two variables at once without noticing
self.p_old = self.p.copy() self.p_old = self.p.copy()
self.v_old = self.v.copy() self.v_old = self.v.copy()

237
Druckrohrleitung/Druckrohrleitung_test.ipynb
View File

@@ -1,237 +0,0 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from Druckrohrleitung_class_file import Druckrohrleitung_class\n",
"import matplotlib.pyplot as plt\n",
"\n",
"#importing pressure conversion function\n",
"import sys\n",
"import os\n",
"current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n",
"parent = os.path.dirname(current)\n",
"sys.path.append(parent)\n",
"from functions.pressure_conversion import pressure_conversion"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib qt5\n",
"#define constants\n",
"\n",
"g = 9.81 # gravitational acceleration [m/s²]\n",
"rho = 1000 # density of water [kg/m³]\n",
"\n",
"L = 1000 # length of pipeline [m]\n",
"D = 1 # pipe diameter [m]\n",
"Q0 = 2 # initial flow in whole pipe [m³/s]\n",
"h_res = 20 # water level in upstream reservoir [m]\n",
"n = 10 # number of pipe segments in discretization\n",
"nt = 100 # number of time steps after initial conditions\n",
"f_D = 0.01 # Darcy friction factor\n",
"c = 400 # propagation velocity of the pressure wave [m/s]\n",
"h_pipe = 200 # hydraulic head without reservoir [m] \n",
"alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n",
"\n",
"\n",
"# preparing the discretization and initial conditions\n",
"\n",
"dx = L/n # length of each pipe segment\n",
"dt = dx/c # timestep according to method of characterisitics\n",
"nn = n+1 # number of nodes\n",
"pl_vec = np.arange(0,nn*dx,dx) # pl = pipe-length. position of the nodes on the pipeline\n",
"t_vec = np.arange(0,nt*dt,dt) # time vector\n",
"h_vec = np.arange(0,h_pipe+h_pipe/n,h_pipe/n) # hydraulic head of pipeline at each node\n",
"\n",
"v_init = np.full(nn,Q0/(D**2/4*np.pi))\n",
"p_init = (rho*g*(h_res+h_vec)-v_init**2*rho/2)-(f_D*pl_vec/D*rho/2*v_init**2) # ref Wikipedia: Darcy Weisbach\n",
"\n",
"# storage vectors for old parameters\n",
"v_old = v_init.copy()\n",
"p_old = p_init.copy() \n",
"\n",
"# storage vectors for new parameters\n",
"v_new = np.empty_like(v_old)\n",
"p_new = np.empty_like(p_old)\n",
"\n",
"# storage vector for time evolution of parameters at node 0 (at reservoir)\n",
"p_0 = np.full_like(t_vec,p_init[0])\n",
"v_0 = np.full_like(t_vec,v_init[0])\n",
"\n",
"# storage vector for time evolution of parameters at node N+1 (at valve)\n",
"p_np1 = np.full_like(t_vec,p_init[-1])\n",
"v_np1 = np.full_like(t_vec,v_init[-1])\n",
"\n",
"for it in range(1,nt):\n",
"\n",
" # set boundary conditions\n",
" v_new[-1] = 0 # in front of the instantaneously closing valve, the velocity is 0\n",
" p_new[0] = p_init[0] # hydrostatic pressure from the reservoir\n",
"\n",
" # calculate the new parameters at first and last node\n",
" v_new[0] = v_old[1]+1/(rho*c)*(p_init[0]-p_old[1])+dt*g*np.sin(alpha)-f_D*dt/(2*D)*abs(v_old[1])*v_old[1]\n",
" p_new[-1] = p_old[-2]+rho*c*v_old[-2]-rho*c*f_D*dt/(2*D) *abs(v_old[-2])*v_old[-2]\n",
"\n",
" # calculate parameters at second to second-to-last nodes \n",
" #equation 2-30 plus 2-31 (and refactor for v_i^j+1) in block 2\n",
"\n",
" for i in range(1,nn-1):\n",
" v_new[i] = 0.5*(v_old[i-1]+v_old[i+1])+0.5/(rho*c)*(p_old[i-1]-p_old[i+1]) \\\n",
" +dt*g*np.sin(alpha)-f_D*dt/(4*D)*(abs(v_old[i-1])*v_old[i-1]+abs(v_old[i+1])*v_old[i+1])\n",
"\n",
" p_new[i] = 0.5*rho*c*(v_old[i-1]-v_old[i+1])+0.5*(p_old[i-1]+p_old[i+1]) \\\n",
" -rho*c*f_D*dt/(4*D)*(abs(v_old[i-1])*v_old[i-1]-abs(v_old[i+1])*v_old[i+1])\n",
" \n",
"\n",
" # prepare for next loop\n",
" # use .copy() to avoid that memory address is overwritten and hell breaks loose :D\n",
" #https://www.geeksforgeeks.org/array-copying-in-python/\n",
" p_old = p_new.copy()\n",
" v_old = v_new.copy()\n",
"\n",
" # store parameters of node 1 (at reservoir)\n",
" p_0[it] = p_new[0]\n",
" v_0[it] = v_new[0]\n",
" # store parameters of node N+1 (at reservoir)\n",
" p_np1[it] = p_new[-1]\n",
" v_np1[it] = v_new[-1]"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n",
"\n",
"pipe.set_pressure_propagation_velocity(c)\n",
"pipe.set_number_of_timesteps(nt)\n",
"\n",
"pipe.set_initial_pressure(p_init)\n",
"pipe.set_initial_flow_velocity(v_init)\n",
"pipe.set_boundary_conditions_next_timestep(v_0[0],p_0[0],v_np1[0])\n",
"\n",
"# storage vector for time evolution of parameters at node 0 (at reservoir)\n",
"pipe.p_0 = np.full_like(t_vec,p_init[0])\n",
"pipe.v_0 = np.full_like(t_vec,v_init[0])\n",
"\n",
"# storage vector for time evolution of parameters at node N+1 (at valve)\n",
"pipe.p_np1 = np.full_like(t_vec,p_init[-1])\n",
"pipe.v_np1 = np.full_like(t_vec,v_init[-1])\n",
"\n",
"fig2,axs2 = plt.subplots(2,1)\n",
"axs2[0].set_title('Pressure distribution in pipeline')\n",
"axs2[1].set_title('Velocity distribution in pipeline')\n",
"axs2[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs2[0].set_ylabel(r'$p$ [mWS]')\n",
"axs2[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs2[1].set_ylabel(r'$p$ [mWS]')\n",
"lo_00, = axs2[0].plot(pl_vec,pressure_conversion(pipe.p_old,'Pa','mWS')[0],marker='.')\n",
"lo_01, = axs2[1].plot(pl_vec,pipe.v_old,marker='.')\n",
"axs2[0].set_ylim([-2*np.max(pressure_conversion(p_init,'Pa','mWS')[0]),2*np.max(pressure_conversion(p_init,'Pa','mWS')[0])])\n",
"axs2[1].set_ylim([-2*np.max(v_init),2*np.max(v_init)])\n",
"fig2.tight_layout()\n",
"\n",
"\n",
"for it in range(1,pipe.nt):\n",
" pipe.set_boundary_conditions_next_timestep(v_0[it],p_0[it],v_np1[it])\n",
" pipe.timestep_characteristic_method()\n",
" lo_00.set_ydata(pressure_conversion(pipe.p,'Pa','mWS')[0])\n",
" lo_01.set_ydata(pipe.v)\n",
"\n",
" # store parameters of node 0 (at reservoir)\n",
" pipe.p_0[it] = pipe.p[0]\n",
" pipe.v_0[it] = pipe.v[0]\n",
" # store parameters of node N+1 (at reservoir)\n",
" pipe.p_np1[it] = pipe.p[-1]\n",
" pipe.v_np1[it] = pipe.v[-1]\n",
" \n",
" fig2.suptitle(str(it))\n",
" fig2.canvas.draw()\n",
" fig2.tight_layout()\n",
" plt.pause(0.2)\n"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"fig3,axs3 = plt.subplots(2,2)\n",
"axs3[0,0].plot(t_vec,pressure_conversion(pipe.p_0,'Pa','mWS')[0])\n",
"axs3[0,1].plot(t_vec,pipe.v_0)\n",
"axs3[1,0].plot(t_vec,pressure_conversion(pipe.p_np1,'Pa','mWS')[0])\n",
"axs3[1,1].plot(t_vec,pipe.v_np1)\n",
"axs3[0,0].set_title('Pressure Reservoir')\n",
"axs3[0,1].set_title('Velocity Reservoir')\n",
"axs3[1,0].set_title('Pressure Turbine')\n",
"axs3[1,1].set_title('Velocity Turbine')\n",
"axs3[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs3[0,0].set_ylabel(r'$p$ [mWS]')\n",
"axs3[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs3[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
"axs3[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs3[1,0].set_ylabel(r'$p$ [mWS]')\n",
"axs3[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs3[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
"fig3.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.29590621205048523\n"
]
}
],
"source": [
"print(np.mean(v_0))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.8.13 ('Georg_DT_Slot3')",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.13"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "84fb123bdc47ab647d3782661abcbe80fbb79236dd2f8adf4cef30e8755eb2cd"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}

202
Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb Normal file
View File

@@ -0,0 +1,202 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from Druckrohrleitung_class_file import Druckrohrleitung_class\n",
"import matplotlib.pyplot as plt\n",
"\n",
"#importing pressure conversion function\n",
"import sys\n",
"import os\n",
"current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n",
"parent = os.path.dirname(current)\n",
"sys.path.append(parent)\n",
"from functions.pressure_conversion import pressure_conversion"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib qt5\n",
"#define constants\n",
"\n",
"g = 9.81 # gravitational acceleration [m/s²]\n",
"rho = 1000 # density of water [kg/m³]\n",
"\n",
"L = 1000 # length of pipeline [m]\n",
"D = 1 # pipe diameter [m]\n",
"Q0 = 2 # initial flow in whole pipe [m³/s]\n",
"h_res = 20 # water level in upstream reservoir [m]\n",
"n = 10 # number of pipe segments in discretization\n",
"nt = 100 # number of time steps after initial conditions\n",
"f_D = 0.01 # Darcy friction factor\n",
"c = 400 # propagation velocity of the pressure wave [m/s]\n",
"h_pipe = 200 # hydraulic head without reservoir [m] \n",
"alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n",
"\n",
"\n",
"# preparing the discretization and initial conditions\n",
"initial_influx = 2. # m³/s\n",
"initial_level = 10. # 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"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"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)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# initialization for timeloop\n",
"\n",
"# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n",
"v_old = pipe.get_current_velocity_distribution()\n",
"p_old = pipe.get_current_pressure_distribution()\n",
"\n",
"# prepare the vectors in which the temporal evolution of the boundary conditions are stored\n",
" # keep in mind, that the velocity at the turbine and the pressure at the reservoir are set manually and\n",
" # through the time evolution of the reservoir respectively \n",
" # the pressure at the turbine and the velocity at the reservoir are calculated from the method of characteristics\n",
"v_boundary_res = np.zeros_like(t_vec)\n",
"v_boundary_tur = np.zeros_like(t_vec)\n",
"p_boundary_res = np.zeros_like(t_vec)\n",
"p_boundary_tur = np.zeros_like(t_vec)\n",
"\n",
"# set the boundary conditions for the first timestep\n",
"v_boundary_res[0] = v_old[0]\n",
"v_boundary_tur[0] = v_old[-1] \n",
"p_boundary_res[0] = p_old[0]\n",
"p_boundary_tur[0] = p_old[-1]\n"
]
},
{
"cell_type": "code",
"execution_count": 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",
"\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",
"\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",
"\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",
"\n",
"fig3.tight_layout()\n",
"plt.show()"
]
}
],
"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
}

0
Regler/regler_test_optimierung.ipynb → Regler/old/regler_test_optimierung.ipynb
View File

8
Turbinen/Turbinen_class_file.py
View File

@@ -12,14 +12,16 @@ class Francis_Turbine:
self.Q_n = Q_nenn self.Q_n = Q_nenn
self.p_n = p_nenn self.p_n = p_nenn
self.LA_n = 1. # 100% self.LA_n = 1. # 100%
h,_ = pressure_conversion(p_nenn,'Pa','MWs') h = pressure_conversion(p_nenn,'Pa','MWs')
self.A = Q_nenn/(np.sqrt(2*9.81*h)*0.98) self.A = Q_nenn/(np.sqrt(2*9.81*h)*0.98)
def set_LA(self,LA): def set_LA(self,LA):
self.LA = LA self.LA = LA
def set_pressure(self,pressure):
self.p = pressure
def get_Q(self,p): def get_Q(self):
self.Q = self.Q_n*(self.LA/self.LA_n)*np.sqrt(p/self.p_n) self.Q = self.Q_n*(self.LA/self.LA_n)*np.sqrt(self.p/self.p_n)
return self.Q return self.Q
def set_closing_time(self,t_closing): def set_closing_time(self,t_closing):

167
Turbinen/Turbinen_test.ipynb
View File

File diff suppressed because one or more lines are too long

40
Untertweng.ipynb
View File

@@ -25,7 +25,7 @@
"\n", "\n",
"#Turbine\n", "#Turbine\n",
"Q_nenn = 0.85\n", "Q_nenn = 0.85\n",
"p_nenn,_ = pressure_conversion(10.6,'bar','Pa')\n", "p_nenn = pressure_conversion(10.6,'bar','Pa')\n",
"\n", "\n",
"# physics\n", "# physics\n",
"g = 9.81 # gravitational acceleration [m/s²]\n", "g = 9.81 # gravitational acceleration [m/s²]\n",
@@ -104,12 +104,12 @@
"# 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_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", "\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_steady_state(initial_influx,V.level,pl_vec,h_vec,initial_pressure_unit,conversion_pressure_unit)\n", "pipe.set_steady_state(initial_influx,V.level,pl_vec,h_vec)\n",
"\n", "\n",
"\n", "\n",
"T1 = Francis_Turbine(Q_nenn,p_nenn)\n", "T1 = Francis_Turbine(Q_nenn,p_nenn)\n",
@@ -137,14 +137,14 @@
" # 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",
" # through the time evolution of the reservoir respectively \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", " # 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_res = np.zeros_like(t_vec)\n",
"v_boundary_tur = np.empty_like(t_vec)\n", "v_boundary_tur = np.zeros_like(t_vec)\n",
"p_boundary_res = np.empty_like(t_vec)\n", "p_boundary_res = np.zeros_like(t_vec)\n",
"p_boundary_tur = np.empty_like(t_vec)\n", "p_boundary_tur = np.zeros_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,V.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.zeros([nt_eRK4]) # level throughout each reservoir timestep-used for plotting and overwritten afterwards\n",
"\n", "\n",
"# set the boundary 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",
@@ -176,7 +176,7 @@
"axs1[0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", "axs1[0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
"axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs1[1].set_ylabel(r'$v$ [$\\mathrm{m} / \\mathrm{s}$]')\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_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_01, = axs1[1].plot(pl_vec,pipe.v_old,marker='.')\n",
"axs1[0].autoscale()\n", "axs1[0].autoscale()\n",
"axs1[1].autoscale()\n", "axs1[1].autoscale()\n",
@@ -193,7 +193,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 6, "execution_count": 7,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -208,7 +208,7 @@
" 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",
" V.level = V.update_level(V.timestep) # \n", " V.level = V.update_level(V.timestep) # \n",
" V.set_volume() # update volume in reservoir\n", " V.update_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", " i_max = it_pipe # for plotting only calculated values\n",
@@ -217,15 +217,14 @@
"\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.outflux_vel**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",
" +dt*g*np.sin(alpha)\n",
"\n", "\n",
" T1.change_LA(LA_soll_vec[it_pipe],dt)\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", " v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_Q(p_old[-1])\n",
"\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(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",
" v_boundary_res[it_pipe] = pipe.v_boundary_res\n",
"\n", "\n",
" # perform the next timestep via the characteristic method\n", " # perform the next timestep via the characteristic method\n",
" pipe.timestep_characteristic_method()\n", " pipe.timestep_characteristic_method()\n",
@@ -236,7 +235,7 @@
" lo_01.remove()\n", " lo_01.remove()\n",
" # lo_02.remove()\n", " # lo_02.remove()\n",
" # plot new pressure and velocity distribution in the pipeline\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_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_01, = axs1[1].plot(pl_vec,pipe.v_old,marker='.',c='blue')\n",
" # lo_02, = axs1[2].plot(level_vec_2,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.suptitle(str(round(t_vec[it_pipe],2))+ ' s / '+str(round(t_vec[-1],2)) + ' s' )\n",
@@ -255,16 +254,17 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 7, "execution_count": 9,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# plot time evolution of boundary pressure and velocity as well as the reservoir level\n", "# plot time evolution of boundary pressure and velocity as well as the reservoir level\n",
"\n", "\n",
"fig2,axs2 = plt.subplots(3,2)\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,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].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[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[1,1].plot(t_vec,v_boundary_tur)\n",
"axs2[2,0].plot(t_vec,level_vec)\n", "axs2[2,0].plot(t_vec,level_vec)\n",
"axs2[0,0].set_title('Pressure reservoir')\n", "axs2[0,0].set_title('Pressure reservoir')\n",
@@ -290,7 +290,7 @@
], ],
"metadata": { "metadata": {
"kernelspec": { "kernelspec": {
"display_name": "Python 3.8.13 ('DT_Slot_3')", "display_name": "Python 3.8.13 ('Georg_DT_Slot3')",
"language": "python", "language": "python",
"name": "python3" "name": "python3"
}, },
@@ -309,7 +309,7 @@
"orig_nbformat": 4, "orig_nbformat": 4,
"vscode": { "vscode": {
"interpreter": { "interpreter": {
"hash": "4a28055eb8a3160fa4c7e4fca69770c4e0a1add985300856aa3fcf4ce32a2c48" "hash": "84fb123bdc47ab647d3782661abcbe80fbb79236dd2f8adf4cef30e8755eb2cd"
} }
} }
}, },