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

Merge branch 'Dev'

Browse Source
This commit is contained in:
Brantegger Georg
2022-08-08 07:53:06 +02:00
parent 1aee8debc4 5a790d5ca5
commit 112c9d427a
11 changed files with 1432 additions and 856 deletions
Download Patch File Download Diff File Expand all files Collapse all files

181
Ausgleichsbecken/Ausgleichsbecken_class_file.py
View File

@@ -17,8 +17,8 @@ def FODE_function(x_out,h,A,A_a,p,rho,g):
# and flow direction # and flow direction
# x_out ... effusion velocity # x_out ... effusion velocity
# h ... level in the reservoir # h ... level in the reservoir
# A_a ... Outflux_Area # A_a ... Area_outflux
# A ... Reservoir_Area # A ... Area_reservoir_base
# g ... gravitational acceleration # g ... gravitational acceleration
# rho ... density of the liquid in the reservoir # rho ... density of the liquid in the reservoir
f = x_out*abs(x_out)/h*(A_a/A-1.)+g-p/(rho*h) f = x_out*abs(x_out)/h*(A_a/A-1.)+g-p/(rho*h)
@@ -28,7 +28,7 @@ def FODE_function(x_out,h,A,A_a,p,rho,g):
class Ausgleichsbecken_class: class Ausgleichsbecken_class:
# units # units
# make sure that units and print units are the same # 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() # units are used to label graphs and disp units are used to have a bearable format when using pythons print()
area_unit = r'$\mathrm{m}^2$' area_unit = r'$\mathrm{m}^2$'
area_outflux_unit = r'$\mathrm{m}^2$' area_outflux_unit = r'$\mathrm{m}^2$'
density_unit = r'$\mathrm{kg}/\mathrm{m}^3$' density_unit = r'$\mathrm{kg}/\mathrm{m}^3$'
@@ -39,157 +39,191 @@ class Ausgleichsbecken_class:
velocity_unit = r'$\mathrm{m}/\mathrm{s}$' velocity_unit = r'$\mathrm{m}/\mathrm{s}$'
volume_unit = r'$\mathrm{m}^3$' volume_unit = r'$\mathrm{m}^3$'
area_unit_print = 'm²' area_unit_disp = 'm²'
area_outflux_unit_print = 'm²' area_outflux_unit_disp = 'm²'
density_unit_print = 'kg/m³' density_unit_disp = 'kg/m³'
flux_unit_print = 'm³/s' flux_unit_disp = 'm³/s'
level_unit_print = 'm' level_unit_disp = 'm'
pressure_unit_print = '--' # will be set by .set_pressure() method time_unit_disp = 's'
time_unit_print = 's' velocity_unit_disp = 'm/s'
velocity_unit_print = 'm/s' volume_unit_disp = 'm³'
volume_unit_print = 'm³'
g = 9.81 # m/s² gravitational acceleration g = 9.81 # m/s² gravitational acceleration
# init # init
def __init__(self,area,outflux_area,level_min = 0,level_max = np.inf ,timestep = 1,rho = 1000): def __init__(self,area,area_outflux,timestep,pressure_unit_disp,level_min=0,level_max=np.inf,rho = 1000.):
self.area = area # base area of the rectangular structure self.area = area # base area of the cuboid reservoir
self.area_outflux = outflux_area # area of the outlet towards the pipeline self.area_out = area_outflux # area of the outlet towards the pipeline
self.density = rho # density of the liquid in the system self.density = rho # density of the liquid in the system
self.level_min = level_min # lowest allowed water level self.level_min = level_min # lowest allowed water level
self.level_max = level_max # highest allowed water level self.level_max = level_max # highest allowed water level
self.timestep = timestep # timestep of the simulation self.pressure_unit_disp = pressure_unit_disp # pressure unit for displaying
self.timestep = timestep # timestep in the time evolution method
# initialize for get_info # initialize for get_info
self.influx = "--" self.influx = "--"
self.level = "--"
self.outflux = "--" self.outflux = "--"
self.level = "--"
self.pressure = "--"
self.volume = "--" self.volume = "--"
# setter # setter
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 initial 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(set_flag=True)
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,set_flag=True) method to update level based on net flux.')
def set_level(self,level): def set_initial_pressure(self,initial_pressure):
self.level = level # sets the initial static pressure present at the outlet of the reservoir and should only be called during initialization
if self.pressure == '--':
self.pressure = initial_pressure
else:
raise Exception('Initial pressure was already set once. Use the .update_pressure(self) method to update pressure based current level.')
def set_influx(self,influx): def set_influx(self,influx):
# sets influx to the reservoir in m³/s # sets influx to the reservoir in m³/s
# positive influx means that liquid flows into the reservoir # positive influx means that liquid flows into the reservoir
self.influx = influx self.influx = influx
def set_outflux(self,outflux): def set_outflux(self,outflux,display_warning=True):
# sets outflux to the reservoir in m³/s # sets outflux to the reservoir in m³/s
# positive outflux means that liquid flows out of reservoir the reservoir # positive outflux means that liquid flows out of reservoir the reservoir
if display_warning == True:
print('You are setting the outflux from the reservoir manually. \n \
This is not an intended use of this method. \n \
Refer to the timestep_reservoir_evolution() method instead.')
self.outflux = outflux self.outflux = outflux
def set_initial_pressure(self,pressure,display_pressure_unit): def set_level(self,level,display_warning=True):
# sets the static pressure present at the outlet of the reservoir # sets level in the reservoir in m
# units are used to convert and display the pressure if display_warning == True:
self.pressure = pressure print('You are setting the level of the reservoir manually. \n \
self.pressure_unit_print = display_pressure_unit This is not an intended use of this method. \n \
Refer to the update_level() method instead.')
self.level = level
def set_pressure(self,pressure): def set_pressure(self,pressure,display_warning=True):
# sets the static pressure present at the outlet of the reservoir # sets pressure in the pipeline just below the reservoir in Pa
if display_warning == True:
print('You are setting the pressure below the reservoir manually. \n \
This is not an intended use of this method. \n \
Refer to the update_pressure() method instead.')
self.pressure = pressure self.pressure = pressure
def set_steady_state(self,ss_influx,ss_level,display_pressure_unit): def set_volume(self,volume,display_warning=True):
if display_warning == True:
print('You are setting the volume in the reservoir manually. \n \
This is not an intended use of this method. \n \
Refer to the .update_volume() or set_initial_level() method instead.')
self.volume = volume
def set_steady_state(self,ss_influx,ss_level):
# set 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 area 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_influx_vel = abs(ss_influx/self.area) ss_influx_vel = abs(ss_influx/self.area)
ss_outflux_vel = abs(ss_outflux/self.area_outflux) ss_outflux_vel = abs(ss_outflux/self.area_out)
ss_pressure = self.density*self.g*ss_level+self.density*ss_outflux_vel*(ss_influx_vel-ss_outflux_vel) ss_pressure = self.density*self.g*ss_level+self.density*ss_outflux_vel*(ss_influx_vel-ss_outflux_vel)
self.set_influx(ss_influx) self.set_influx(ss_influx)
self.set_initial_level(ss_level) self.set_initial_level(ss_level)
self.set_initial_pressure(ss_pressure,display_pressure_unit) self.set_initial_pressure(ss_pressure)
self.set_outflux(ss_outflux) self.set_outflux(ss_outflux,display_warning=False)
# 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_disp)
outflux_vel = self.outflux/self.area_outflux outflux_vel = self.outflux/self.area_out
if full == True: if full == True:
# :<10 pads the self.value to be 10 characters wide # :<10 pads the self.value to be 10 characters wide
print_str = (f"The cuboid reservoir has the following attributes: {new_line}" print_str = (f"The cuboid reservoir has the following attributes: {new_line}"
f"----------------------------- {new_line}" f"----------------------------- {new_line}"
f"Base area = {self.area:<10} {self.area_unit_print} {new_line}" f"Base area = {self.area:<10} {self.area_unit_disp} {new_line}"
f"Outflux area = {round(self.area_outflux,3):<10} {self.area_outflux_unit_print} {new_line}" f"Outflux area = {round(self.area_out,3):<10} {self.area_out_unit_disp} {new_line}"
f"Current level = {self.level:<10} {self.level_unit_print}{new_line}" f"Current level = {self.level:<10} {self.level_unit_disp}{new_line}"
f"Critical level low = {self.level_min:<10} {self.level_unit_print} {new_line}" f"Critical level low = {self.level_min:<10} {self.level_unit_disp} {new_line}"
f"Critical level high = {self.level_max:<10} {self.level_unit_print} {new_line}" f"Critical level high = {self.level_max:<10} {self.level_unit_disp} {new_line}"
f"Volume in reservoir = {self.volume:<10} {self.volume_unit_print} {new_line}" f"Volume in reservoir = {self.volume:<10} {self.volume_unit_disp} {new_line}"
f"Current influx = {self.influx:<10} {self.flux_unit_print} {new_line}" f"Current influx = {self.influx:<10} {self.flux_unit_disp} {new_line}"
f"Current outflux = {self.outflux:<10} {self.flux_unit_print} {new_line}" f"Current outflux = {self.outflux:<10} {self.flux_unit_disp} {new_line}"
f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_print} {new_line}" f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_disp} {new_line}"
f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}" f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_disp} {new_line}"
f"Simulation timestep = {self.timestep:<10} {self.time_unit_print} {new_line}" f"Simulation timestep = {self.timestep:<10} {self.time_unit_disp} {new_line}"
f"Density of liquid = {self.density:<10} {self.density_unit_print} {new_line}" f"Density of liquid = {self.density:<10} {self.density_unit_disp} {new_line}"
f"----------------------------- {new_line}") f"----------------------------- {new_line}")
else: else:
# :<10 pads the self.value to be 10 characters wide # :<10 pads the self.value to be 10 characters wide
print_str = (f"The current attributes are: {new_line}" print_str = (f"The current attributes are: {new_line}"
f"----------------------------- {new_line}" f"----------------------------- {new_line}"
f"Current level = {self.level:<10} {self.level_unit_print}{new_line}" f"Current level = {self.level:<10} {self.level_unit_disp}{new_line}"
f"Volume in reservoir = {self.volume:<10} {self.volume_unit_print} {new_line}" f"Current volume = {self.volume:<10} {self.volume_unit_disp} {new_line}"
f"Current influx = {self.influx:<10} {self.flux_unit_print} {new_line}" f"Current influx = {self.influx:<10} {self.flux_unit_disp} {new_line}"
f"Current outflux = {self.outflux:<10} {self.flux_unit_print} {new_line}" f"Current outflux = {self.outflux:<10} {self.flux_unit_disp} {new_line}"
f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_print} {new_line}" f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_disp} {new_line}"
f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}" f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_disp} {new_line}"
f"----------------------------- {new_line}") f"----------------------------- {new_line}")
print(print_str) print(print_str)
def get_current_level(self):
return self.level
def get_current_influx(self): def get_current_influx(self):
return self.influx return self.influx
def get_current_outflux(self): def get_current_outflux(self):
return self.outflux return self.outflux
def get_current_level(self):
return self.level
def get_current_pressure(self): def get_current_pressure(self):
return self.pressure return self.pressure
def get_current_volume(self):
return self.volume
# update methods
# methods def update_level(self,timestep,set_flag=False):
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
# 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 net_flux = self.influx-self.outflux
delta_level = net_flux*timestep/self.area delta_level = net_flux*timestep/self.area
new_level = (self.level+delta_level) level_new = (self.level+delta_level)
return new_level if set_flag == True:
self.set_level(level_new,display_warning=False)
elif set_flag == False:
return level_new
def update_pressure(self): def update_pressure(self,set_flag=False):
influx_vel = abs(self.influx/self.area) influx_vel = abs(self.influx/self.area)
outflux_vel = abs(self.outflux/self.area_outflux) outflux_vel = abs(self.outflux/self.area_out)
p_new = self.density*self.g*self.level+self.density*outflux_vel*(influx_vel-outflux_vel) p_new = self.density*self.g*self.level+self.density*outflux_vel*(influx_vel-outflux_vel)
if set_flag ==True:
self.set_pressure(p_new,display_warning=False)
elif set_flag == False:
return p_new return p_new
def update_volume(self,set_flag=False):
volume_new = self.level*self.area
if set_flag == True:
self.set_volume(volume_new,display_warning=False)
elif set_flag == False:
return volume_new
#methods
def timestep_reservoir_evolution(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
dt = self.timestep dt = self.timestep
rho = self.density rho = self.density
g = self.g g = self.g
A = self.area A = self.area
A_a = self.area_outflux A_a = self.area_out
yn = self.outflux/A_a # outflux velocity yn = self.outflux/A_a # outflux velocity
h = self.level h = self.level
h_hs = self.update_level(dt/2) h_hs = self.update_level(dt/2)
@@ -203,10 +237,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))
new_outflux = ynp1*A_a self.set_outflux(ynp1*A_a,display_warning=False)
new_level = self.update_level(dt) self.update_level(dt,set_flag=True)
new_pressure = self.update_pressure() self.update_volume(set_flag=True)
self.update_pressure(set_flag=True)
self.set_outflux(new_outflux)
self.set_level(new_level)
self.set_pressure(new_pressure)

187
Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 1, "execution_count": 29,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -21,126 +21,137 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 2, "execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"L = 1000.\n",
"n = 10000 # number of pipe segments in discretization\n",
"c = 400. \n",
"dx = L/n # length of each pipe segment\n",
"dt = dx/c \n",
"\n",
"# # define constants\n",
"# initial_level = 10.1 # m\n",
"# initial_influx = 0.8 # m³/s\n",
"# conversion_pressure_unit = 'mWS'\n",
"\n",
"# area_base = 75. # m²\n",
"# area_outflux = (0.9/2)**2*np.pi # m²\n",
"# critical_level_low = 0. # m\n",
"# critical_level_high = 10. # m\n",
"# simulation_timestep = dt # s\n",
"\n",
"# # for while loop\n",
"# total_min_level = 0.01 # m\n",
"# total_max_time = 100 # s\n",
"\n",
"# nt = int(total_max_time//simulation_timestep)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# define constants\n", "# define constants\n",
"initial_level = 10.1 # m\n",
"initial_influx = 1. # m³/s\n",
"# initial_outflux = 1. # m³/s\n",
"# initial_pipeline_pressure = 10.\n",
"# initial_pressure_unit = 'mWS'\n",
"conversion_pressure_unit = 'mWS'\n",
"\n", "\n",
"area_base = 75. # m²\n", " # for physics\n",
"area_outflux = 2. # m²\n", "g = 9.81 # [m/s²] gravitational acceleration \n",
"critical_level_low = 0. # m\n", "rho = 1000. # [kg/m³] density of water \n",
"critical_level_high = 10. # m\n", "pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n",
"simulation_timestep = dt # s\n", "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n",
"\n", "\n",
"# for while loop\n",
"total_min_level = 0.01 # m\n",
"total_max_time = 100 # s\n",
"\n", "\n",
"nt = int(total_max_time//simulation_timestep)" " # for Turbine\n",
"Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n",
"Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n",
"Tur_closingTime = 90. # [s] closing time of turbine\n",
"\n",
"\n",
" # for PI controller\n",
"Con_targetLevel = 8. # [m]\n",
"Con_K_p = 0.1 # [-] proportional constant of PI controller\n",
"Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\n",
"Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n",
"\n",
"\n",
" # for pipeline\n",
"Pip_length = (535.+478.) # [m] length of pipeline\n",
"Pip_dia = 0.9 # [m] diameter of pipeline\n",
"Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n",
"Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n",
"Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n",
"Pip_n_seg = 50 # [-] number of pipe segments in discretization\n",
"Pip_f_D = 0.014 # [-] Darcy friction factor\n",
"Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n",
" # derivatives of the pipeline constants\n",
"Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n",
"Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n",
"Pip_nn = Pip_n_seg+1 # [1] number of nodes\n",
"Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n",
"Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n",
"\n",
"\n",
" # for reservoir\n",
"Res_area_base = 5. # [m²] total base are of the cuboid reservoir \n",
"Res_area_out = Pip_area # [m²] outflux area of the reservoir, given by pipeline area\n",
"Res_level_crit_lo = 0. # [m] for yet-to-be-implemented warnings\n",
"Res_level_crit_hi = np.inf # [m] for yet-to-be-implemented warnings\n",
"Res_dt_approx = 1e-3 # [s] approx. timestep of reservoir time evolution to ensure numerical stability (see Res_nt why approx.)\n",
"Res_nt = max(1,int(Pip_dt//Res_dt_approx)) # [1] number of timesteps of the reservoir time evolution within one timestep of the pipeline\n",
"Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n",
"\n",
" # for general simulation\n",
"flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n",
"level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n",
"simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n",
"nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n",
"t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 4, "execution_count": 31,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"%matplotlib qt\n", "# create objects\n",
"\n", "\n",
"V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,simulation_timestep)\n", "# Upstream reservoir\n",
"# V.set_initial_level(initial_level) \n", "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n",
"# V.set_influx(initial_influx)\n", "reservoir.set_steady_state(flux_init,level_init)\n",
"# V.set_outflux(initial_outflux)\n",
"# V.set_initial_pressure(pressure_conversion(initial_pipeline_pressure,input_unit = initial_pressure_unit, target_unit = 'Pa'),conversion_pressure_unit)\n",
"# V.pressure = converted_pressure\n",
"V.set_steady_state(initial_influx,initial_level,conversion_pressure_unit)\n",
"\n", "\n",
"time_vec = np.arange(0,nt+1,1)*simulation_timestep\n", "reservoir.get_info(full=True)\n",
"outflux_vec = np.zeros_like(time_vec)\n",
"outflux_vec[0] = V.get_current_outflux()\n",
"level_vec = np.zeros_like(time_vec)\n",
"level_vec[0] = V.get_current_level()\n",
"pressure_vec = np.zeros_like(time_vec)\n",
"pressure_vec[0] = V.get_current_pressure()\n",
"\n",
"# pressure_vec = np.full_like(time_vec,converted_pressure)*((np.sin(time_vec)+1)*np.exp(-time_vec/50))\n",
" \n",
"i_max = -1\n",
"\n", "\n",
"# initialize vectors\n",
"outflux_vec = np.zeros_like(t_vec)\n",
"outflux_vec[0] = reservoir.get_current_outflux()\n",
"level_vec = np.zeros_like(t_vec)\n",
"level_vec[0] = reservoir.get_current_level()\n",
"volume_vec = np.zeros_like(t_vec)\n",
"volume_vec[0] = reservoir.get_current_volume()\n",
"pressure_vec = np.zeros_like(t_vec)\n",
"pressure_vec[0] = reservoir.get_current_pressure()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# time loop\n",
"for i in range(1,nt+1):\n", "for i in range(1,nt+1):\n",
" V.set_pressure(pressure_vec[i-1])\n", " # if i == 500:\n",
" V.set_outflux(outflux_vec[i-1])\n", " # reservoir.set_influx(0.)\n",
" V.timestep_reservoir_evolution()\n", " reservoir.set_pressure(pressure_vec[i-1],display_warning=False)\n",
" outflux_vec[i] = V.get_current_outflux()\n", " reservoir.set_outflux(outflux_vec[i-1],display_warning=False)\n",
" level_vec[i] = V.get_current_level()\n", " for it_res in range(Res_nt):\n",
" pressure_vec[i] = V.get_current_pressure()\n", " reservoir.timestep_reservoir_evolution() \n",
" if V.level < total_min_level:\n", " \n",
" i_max = i\n", " outflux_vec[i] = reservoir.get_current_outflux()\n",
" break\n", " level_vec[i] = reservoir.get_current_level()\n",
"\n" " pressure_vec[i] = reservoir.get_current_pressure()\n",
"\n",
" reservoir.get_info()"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 5, "execution_count": 32,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"\n", "%matplotlib qt5\n",
"fig1, (ax1, ax2, ax3) = plt.subplots(3, 1)\n", "fig1, (ax1, ax2, ax3) = plt.subplots(3, 1)\n",
"fig1.set_figheight(10)\n", "fig1.set_figheight(10)\n",
"fig1.suptitle('Ausgleichsbecken')\n", "fig1.suptitle('Ausgleichsbecken')\n",
"\n", "\n",
"ax1.plot(time_vec[:i_max],level_vec[:i_max], label='Water level')\n", "ax1.plot(t_vec,level_vec, label='Water level')\n",
"ax1.set_ylabel(r'$h$ ['+V.level_unit+']')\n", "ax1.set_ylabel(r'$h$ ['+reservoir.level_unit+']')\n",
"ax1.set_xlabel(r'$t$ ['+V.time_unit+']')\n", "ax1.set_xlabel(r'$t$ ['+reservoir.time_unit+']')\n",
"ax1.legend()\n", "ax1.legend()\n",
"\n", "\n",
"ax2.plot(time_vec[:i_max],outflux_vec[:i_max], label='Outflux')\n", "ax2.plot(t_vec,outflux_vec, label='Outflux')\n",
"ax2.set_ylabel(r'$Q_{out}$ ['+V.flux_unit+']')\n", "ax2.set_ylabel(r'$Q_{out}$ ['+reservoir.flux_unit+']')\n",
"ax2.set_xlabel(r'$t$ ['+V.time_unit+']')\n", "ax2.set_xlabel(r'$t$ ['+reservoir.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), label='Pipeline pressure at reservoir')\n", "ax3.plot(t_vec,pressure_conversion(pressure_vec,'Pa',pUnit_conv), label='Pipeline pressure at reservoir')\n",
"ax3.set_ylabel(r'$p_{pipeline}$ ['+conversion_pressure_unit+']')\n", "ax3.set_ylabel(r'$p_{pipeline}$ ['+pUnit_conv+']')\n",
"ax3.set_xlabel(r'$t$ ['+V.time_unit+']')\n", "ax3.set_xlabel(r'$t$ ['+reservoir.time_unit+']')\n",
"ax3.legend()\n", "ax3.legend()\n",
"\n", "\n",
"\n", "\n",

139
Druckrohrleitung/Druckrohrleitung_class_file.py
View File

@@ -1,5 +1,13 @@
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$'
@@ -13,72 +21,66 @@ class Druckrohrleitung_class:
velocity_unit = r'$\mathrm{m}/\mathrm{s}$' # for flux and pressure propagation velocity_unit = r'$\mathrm{m}/\mathrm{s}$' # for flux and pressure propagation
volume_unit = r'$\mathrm{m}^3$' volume_unit = r'$\mathrm{m}^3$'
acceleration_unit_print = 'm/s²' acceleration_unit_disp = 'm/s²'
angle_unit_print = 'rad' angle_unit_disp = 'rad'
area_unit_print = 'm²' area_unit_disp = 'm²'
density_unit_print = 'kg/m³' density_unit_disp = 'kg/m³'
flux_unit_print = 'm³/s' flux_unit_disp = 'm³/s'
length_unit_print = 'm' length_unit_disp = 'm'
time_unit_print = 's' time_unit_disp = 's'
velocity_unit_print = 'm/s' # for flux and pressure propagation velocity_unit_disp = 'm/s' # for flux and pressure propagation
volume_unit_print = 'm³' volume_unit_disp = 'm³'
g = 9.81
# 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,pw_vel,timestep,pressure_unit_disp,rho=1000):
self.length = total_length # total length of the pipeline self.length = total_length # total length of the pipeline
self.dia = diameter # diameter of the pipeline self.dia = diameter # diameter of the pipeline
self.n_seg = number_segments # number of segments for the method of characteristics self.n_seg = number_segments # number of segments for the method of characteristics
self.angle = pipeline_angle # angle of the pipeline 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.c = pw_vel
self.dt = timestep
self.density = rho # density of the liquid in the pipeline self.density = rho # density of the liquid in the pipeline
self.g = g # gravitational acceleration
self.A = (diameter/2)**2*np.pi self.A = (diameter/2)**2*np.pi
self.dx = total_length/number_segments # length of each segment 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 self.x_vec = np.arange(0,(number_segments+1),1)*self.dx # vector giving the distance from each node to the start of the pipeline
# initialize for get_info method self.pressure_unit_disp = pressure_unit_disp
self.c = '--'
self.dt = '--'
# setter # setter
def set_pressure_propagation_velocity(self,c):
self.c = c # propagation velocity of the pressure wave
self.dt = self.dx/c # timestep derived from c, demanded by the method of characteristics
def set_number_of_timesteps(self,number_timesteps):
self.nt = number_timesteps # number of timesteps
if self.c == '--':
raise Exception('Please set the pressure propagation velocity before setting the number of timesteps.')
else:
self.t_vec = np.arange(0,self.nt*self.dt,self.dt)
def set_initial_pressure(self,pressure): def set_initial_pressure(self,pressure):
# initialize the pressure distribution in the pipeline # initialize the pressure distribution in the pipeline
if np.size(pressure) == 1: if np.size(pressure) == 1:
self.p0 = np.full_like(self.l_vec,pressure) p0 = np.full_like(self.x_vec,pressure)
elif np.size(pressure) == np.size(self.l_vec): elif np.size(pressure) == np.size(self.x_vec):
self.p0 = pressure 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.x_vec))
#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 = p0.copy()
self.p = self.p0.copy() self.p = p0.copy()
self.p_min = p0.copy()
self.p_max = p0.copy()
def set_initial_flow_velocity(self,velocity): def set_initial_flow_velocity(self,velocity):
# initialize the velocity distribution in the pipeline # initialize the velocity distribution in the pipeline
if np.size(velocity) == 1: if np.size(velocity) == 1:
self.v0 = np.full_like(self.l_vec,velocity) v0 = np.full_like(self.x_vec,velocity)
elif np.size(velocity) == np.size(self.l_vec): elif np.size(velocity) == np.size(self.x_vec):
self.v0 = velocity v0 = velocity
else: else:
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.x_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 = v0.copy()
self.v = self.v0.copy() self.v = v0.copy()
self.v_min = v0.copy()
self.v_max = v0.copy()
def set_boundary_conditions_next_timestep(self,p_reservoir,v_turbine): def set_boundary_conditions_next_timestep(self,p_reservoir,v_turbine):
# derived from the method of characteristics, one can set the boundary conditions for the pressures and flow velocities at the reservoir and the turbine # derived from the method of characteristics, one can set the boundary conditions for the pressures and flow velocities at the reservoir and the turbine
@@ -112,16 +114,16 @@ class Druckrohrleitung_class:
self.p[0] = p_boundary_res self.p[0] = p_boundary_res
self.p[-1] = p_boundary_tur self.p[-1] = p_boundary_tur
def set_steady_state(self,ss_flux,ss_level_reservoir,area_reservoir,pl_vec,h_vec): def set_steady_state(self,ss_flux,ss_level_reservoir,area_reservoir,x_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 # 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 # the flow velocity is given by the constant flow through the pipe
ss_v0 = np.full(self.n_seg+1,ss_flux/self.A) ss_v0 = np.full_like(self.x_vec,ss_flux/self.A)
# the static pressure is given by static state pressure of the reservoir, corrected for the hydraulic head of the pipe and friction losses # the static pressure is given by static state pressure of the reservoir, corrected for the hydraulic head of the pipe and friction losses
ss_v_in_res = abs(ss_flux/area_reservoir) ss_v_in_res = abs(ss_flux/area_reservoir)
ss_v_out_res = abs(ss_flux/self.A) ss_v_out_res = abs(ss_flux/self.A)
ss_pressure_res = self.density*self.g*(ss_level_reservoir)+self.density*ss_v_out_res*(ss_v_in_res-ss_v_out_res) ss_pressure_res = self.density*self.g*(ss_level_reservoir)+self.density*ss_v_out_res*(ss_v_in_res-ss_v_out_res)
ss_pressure = ss_pressure_res+(self.density*self.g*h_vec)-(self.f_D*pl_vec/self.dia*self.density/2*ss_v0**2) ss_pressure = ss_pressure_res+(self.density*self.g*h_vec)-(self.f_D*x_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) self.set_initial_pressure(ss_pressure)
@@ -135,30 +137,61 @@ class Druckrohrleitung_class:
# :<10 pads the self.value to be 10 characters wide # :<10 pads the self.value to be 10 characters wide
print_str = (f"The pipeline has the following attributes: {new_line}" print_str = (f"The pipeline has the following attributes: {new_line}"
f"----------------------------- {new_line}" f"----------------------------- {new_line}"
f"Length = {self.length:<10} {self.length_unit_print} {new_line}" f"Length = {self.length:<10} {self.length_unit_disp} {new_line}"
f"Diameter = {self.dia:<10} {self.length_unit_print} {new_line}" f"Diameter = {self.dia:<10} {self.length_unit_disp} {new_line}"
f"Number of segments = {self.n_seg:<10} {new_line}" f"Number of segments = {self.n_seg:<10} {new_line}"
f"Number of nodes = {self.n_seg+1:<10} {new_line}" f"Number of nodes = {self.n_seg+1:<10} {new_line}"
f"Length per segments = {self.dx:<10} {self.length_unit_print} {new_line}" f"Length per segments = {self.dx:<10} {self.length_unit_disp} {new_line}"
f"Pipeline angle = {round(self.angle,3):<10} {self.angle_unit_print} {new_line}" f"Pipeline angle = {round(self.angle,3):<10} {self.angle_unit_disp} {new_line}"
f"Pipeline angle = {angle_deg}° {new_line}" f"Pipeline angle = {angle_deg}° {new_line}"
f"Darcy friction factor = {self.f_D:<10} {new_line}" f"Darcy friction factor = {self.f_D:<10} {new_line}"
f"Density of liquid = {self.density:<10} {self.density_unit_print} {new_line}" f"Density of liquid = {self.density:<10} {self.density_unit_disp} {new_line}"
f"Pressure wave vel. = {self.c:<10} {self.velocity_unit_print} {new_line}" f"Pressure wave vel. = {self.c:<10} {self.velocity_unit_disp} {new_line}"
f"Simulation timestep = {self.dt:<10} {self.time_unit_print} {new_line}" f"Simulation timestep = {self.dt:<10} {self.time_unit_disp} {new_line}"
f"Number of timesteps = {self.nt:<10} {new_line}" f"Number of timesteps = {self.nt:<10} {new_line}"
f"Total simulation time = {self.nt*self.dt:<10} {self.time_unit_print} {new_line}" f"Total simulation time = {self.nt*self.dt:<10} {self.time_unit_disp} {new_line}"
f"----------------------------- {new_line}" f"----------------------------- {new_line}"
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_current_pressure_distribution(self): def get_current_pressure_distribution(self,disp=False):
if disp == True:
return pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_disp)
elif disp == False:
return self.p return self.p
def get_current_velocity_distribution(self): def get_current_velocity_distribution(self):
return self.v return self.v
def get_current_flux_distribution(self):
return self.v*self.A
def get_lowest_pressure_per_node(self,disp=False):
if disp == True:
return pressure_conversion(self.p_min,self.pressure_unit,self.pressure_unit_disp)
elif disp == False:
return self.p_min
def get_highest_pressure_per_node(self,disp=False):
if disp == True:
return pressure_conversion(self.p_max,self.pressure_unit,self.pressure_unit_disp)
elif disp == False:
return self.p_max
def get_lowest_velocity_per_node(self):
return self.v_min
def get_highest_velocity_per_node(self):
return self.v_max
def get_lowest_flux_per_node(self):
return self.v_min*self.A
def get_highest_flux_per_node(self):
return self.v_max*self.A
def timestep_characteristic_method(self): def timestep_characteristic_method(self):
# use the method of characteristics to calculate the pressure and velocities at all nodes except the boundary ones # use the method of characteristics to calculate the pressure and velocities at all nodes except the boundary ones
# they are set with the .set_boundary_conditions_next_timestep() method beforehand # they are set with the .set_boundary_conditions_next_timestep() method beforehand
@@ -180,6 +213,12 @@ 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])
# update overall min and max values for pressure and velocity per node
self.p_min = np.minimum(self.p_min,self.p)
self.p_max = np.maximum(self.p_max,self.p)
self.v_min = np.minimum(self.v_min,self.v)
self.v_max = np.maximum(self.v_max,self.v)
# prepare for next call # prepare for next call
# use .copy() to write data to another memory location and avoid the usual python reference pointer # 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 # else one can overwrite data by accidient and change two variables at once without noticing

180
Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb
View File

@@ -26,45 +26,60 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"%matplotlib qt5\n", "# define constants\n",
"#define constants pipe\n",
"\n", "\n",
"g = 9.81 # gravitational acceleration [m/s²]\n", " # for physics\n",
"rho = 1000. # density of water [kg/m³]\n", "g = 9.81 # [m/s²] gravitational acceleration \n",
"\n", "rho = 1000. # [kg/m³] density of water \n",
"L = 1000. # length of pipeline [m]\n", "pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n",
"D = 0.9 # pipe diameter [m]\n", "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n",
"h_res = 10. # water level in upstream reservoir [m]\n",
"n = 50 # number of pipe segments in discretization\n",
"nt = 9000 # number of time steps after initial conditions\n",
"f_D = 0.01 # Darcy friction factor\n",
"c = 400. # propagation velocity of the pressure wave [m/s]\n",
"h_pipe = 105. # hydraulic head without reservoir [m] \n",
"alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n",
"\n", "\n",
"\n", "\n",
"# preparing the discretization and initial conditions\n", " # for Turbine\n",
"initial_flux = 0.8 # m³/s\n", "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n",
"initial_level = h_res # m\n", "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n",
"dx = L/n # length of each pipe segment\n", "Tur_closingTime = 90. # [s] closing time of turbine\n",
"dt = dx/c # timestep according to method of characterisitics\n",
"nn = n+1 # number of nodes\n",
"pl_vec = np.arange(0,nn,1)*dx # pl = pipe-length. position of the nodes on the pipeline\n",
"t_vec = np.arange(0,nt,1)*dt # time vector\n",
"h_vec = np.arange(0,nn,1)*h_pipe/n # hydraulic head of pipeline at each node\n",
"\n", "\n",
"\n", "\n",
"# define constants reservoir\n", " # for PI controller\n",
"conversion_pressure_unit = 'mWS'\n", "Con_targetLevel = 8. # [m]\n",
"Con_K_p = 0.1 # [-] proportional constant of PI controller\n",
"Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\n",
"Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n",
"\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", "\n",
"# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n", " # for pipeline\n",
"nt_eRK4 = 100 # number of simulation steps of reservoir in between timesteps of pipeline \n", "Pip_length = (535.+478.) # [m] length of pipeline\n",
"simulation_timestep = dt/nt_eRK4" "Pip_dia = 0.9 # [m] diameter of pipeline\n",
"Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n",
"Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n",
"Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n",
"Pip_n_seg = 50 # [-] number of pipe segments in discretization\n",
"Pip_f_D = 0.014 # [-] Darcy friction factor\n",
"Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n",
" # derivatives of the pipeline constants\n",
"Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n",
"Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n",
"Pip_nn = Pip_n_seg+1 # [1] number of nodes\n",
"Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n",
"Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n",
"\n",
"\n",
" # for reservoir\n",
"Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n",
"Res_area_out = Pip_area # [m²] outflux area of the reservoir, given by pipeline area\n",
"Res_level_crit_lo = 0. # [m] for yet-to-be-implemented warnings\n",
"Res_level_crit_hi = np.inf # [m] for yet-to-be-implemented warnings\n",
"Res_dt_approx = 1e-3 # [s] approx. timestep of reservoir time evolution to ensure numerical stability (see Res_nt why approx.)\n",
"Res_nt = max(1,int(Pip_dt//Res_dt_approx)) # [1] number of timesteps of the reservoir time evolution within one timestep of the pipeline\n",
"Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n",
"\n",
" # for general simulation\n",
"flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n",
"level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n",
"simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n",
"nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n",
"t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n"
] ]
}, },
{ {
@@ -73,51 +88,37 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"V = Ausgleichsbecken_class(area_base,area_pipe,critical_level_low,critical_level_high,simulation_timestep)\n", "# create objects\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", "# Upstream reservoir\n",
"pipe.set_pressure_propagation_velocity(c)\n", "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n",
"pipe.set_number_of_timesteps(nt)\n", "reservoir.set_steady_state(flux_init,level_init)\n",
"pipe.set_steady_state(initial_flux,initial_level,area_base,pl_vec,h_vec)" "\n",
"# pipeline\n",
"pipe = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_n_seg,Pip_angle,Pip_f_D,Pip_pw_vel,Pip_dt,pUnit_conv,rho)\n",
"pipe.set_steady_state(flux_init,level_init,Res_area_base,Pip_x_vec,Pip_h_vec)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"reservoir.get_info(full=True)\n",
"pipe.get_info(full=True)"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 4, "execution_count": 4,
"metadata": {}, "metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The current attributes are: \n",
"----------------------------- \n",
"Current level = 10.0 m\n",
"Volume in reservoir = -- m³ \n",
"Current influx = 0.8 m³/s \n",
"Current outflux = 0.8 m³/s \n",
"Current outflux vel = 1.258 m/s \n",
"Current pipe pressure = 9.844 mWS \n",
"----------------------------- \n",
"\n"
]
}
],
"source": [
"V.get_info()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# initialization for timeloop\n", "# initialization for timeloop\n",
"\n", "\n",
"level_vec = np.zeros_like(t_vec)\n", "level_vec = np.zeros_like(t_vec)\n",
"level_vec[0] = V.get_current_level()\n", "level_vec[0] = reservoir.get_current_level()\n",
"\n", "\n",
"# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n", "# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n",
"v_old = pipe.get_current_velocity_distribution()\n", "v_old = pipe.get_current_velocity_distribution()\n",
@@ -141,21 +142,23 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 6, "execution_count": 5,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"%matplotlib qt5\n",
"fig1,axs1 = plt.subplots(2,1)\n", "fig1,axs1 = plt.subplots(2,1)\n",
"axs1[0].set_title('Pressure distribution in pipeline')\n", "axs1[0].set_title('Pressure distribution in pipeline')\n",
"axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs1[0].set_ylabel(r'$p$ [mWS]')\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',pUnit_conv)),1.1*np.max(pressure_conversion(p_old,'Pa',pUnit_conv))])\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", "lo_00, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,'Pa',pUnit_conv),marker='.')\n",
"\n", "\n",
"axs1[1].set_title('Velocity distribution in pipeline')\n", "axs1[1].set_title('Velocity distribution in pipeline')\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$ [m/s]')\n", "axs1[1].set_ylabel(r'$v$ [m/s]')\n",
"lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n", "lo_01, = axs1[1].plot(Pip_x_vec,v_old,marker='.')\n",
"axs1[1].autoscale()\n",
"# axs1[1].set_ylim([0.9*np.min(v_old),1.1*np.max(v_boundary_res)])\n", "# axs1[1].set_ylim([0.9*np.min(v_old),1.1*np.max(v_boundary_res)])\n",
"\n", "\n",
"fig1.tight_layout()\n", "fig1.tight_layout()\n",
@@ -164,7 +167,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 7, "execution_count": 6,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -201,21 +204,20 @@
} }
], ],
"source": [ "source": [
"\n", "for it_pipe in range(1,nt+1):\n",
"for it_pipe in range(1,nt):\n",
"# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\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", " # set initial conditions for the reservoir time evolution calculted with e-RK4\n",
" V.set_pressure(p_old[0])\n", " reservoir.set_pressure(p_old[0],display_warning=False)\n",
" V.set_outflux(v_old[0]*area_pipe)\n", " reservoir.set_outflux(v_old[0]*Pip_area,display_warning=False)\n",
" # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n", " # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n",
" for it_res in range(nt_eRK4):\n", " for it_res in range(Res_nt):\n",
" V.timestep_reservoir_evolution() \n", " reservoir.timestep_reservoir_evolution() \n",
" level_vec[it_pipe] = V.get_current_level() \n", " level_vec[it_pipe] = reservoir.get_current_level() \n",
"\n", "\n",
" \n", " \n",
" # set boundary conditions for the next timestep of the characteristic method\n", " # set boundary conditions for the next timestep of the characteristic method\n",
" p_boundary_res[it_pipe] = V.get_current_pressure()\n", " p_boundary_res[it_pipe] = reservoir.get_current_pressure()\n",
" v_boundary_tur[it_pipe] = initial_flux/area_pipe\n", " v_boundary_tur[it_pipe] = flux_init/Pip_area\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(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",
@@ -235,28 +237,30 @@
" 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(p_old,'Pa', conversion_pressure_unit),marker='.',c='blue')\n", " lo_00, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,'Pa', pUnit_conv),marker='.',c='blue')\n",
" lo_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n", " lo_01, = axs1[1].plot(Pip_x_vec,v_old,marker='.',c='blue')\n",
" \n", " \n",
" fig1.suptitle(str(round(t_vec[it_pipe],2)) + '/' + str(round(t_vec[-1],2)))\n", " fig1.suptitle(str(round(t_vec[it_pipe],2)) + '/' + str(round(t_vec[-1],2)))\n",
" fig1.canvas.draw()\n", " fig1.canvas.draw()\n",
" fig1.tight_layout()\n", " fig1.tight_layout()\n",
" plt.pause(0.000001)\n", " plt.pause(0.000001)\n",
"\n" "\n",
"reservoir.get_info(full=True)\n",
"pipe.get_info(full=True)"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 12,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"fig2,axs2 = plt.subplots(2,2)\n", "fig2,axs2 = plt.subplots(2,2)\n",
"axs2[0,0].set_title('Pressure Reservoir')\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].plot(t_vec,pressure_conversion(p_boundary_res,pUnit_calc,pUnit_conv))\n",
"axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[0,0].set_ylabel(r'$p$ [mWS]')\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", "axs2[0,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_res,pUnit_calc,pUnit_conv)),1.1*np.max(pressure_conversion(p_boundary_res,pUnit_calc,pUnit_conv))])\n",
"\n", "\n",
"axs2[0,1].set_title('Velocity Reservoir')\n", "axs2[0,1].set_title('Velocity Reservoir')\n",
"axs2[0,1].plot(t_vec,v_boundary_res)\n", "axs2[0,1].plot(t_vec,v_boundary_res)\n",
@@ -265,16 +269,16 @@
"axs2[0,1].set_ylim([0.9*np.min(v_boundary_res),1.1*np.max(v_boundary_res)])\n", "axs2[0,1].set_ylim([0.9*np.min(v_boundary_res),1.1*np.max(v_boundary_res)])\n",
"\n", "\n",
"axs2[1,0].set_title('Pressure Turbine')\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].plot(t_vec,pressure_conversion(p_boundary_tur,pUnit_calc,pUnit_conv))\n",
"axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[1,0].set_ylabel(r'$p$ [mWS]')\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", "axs2[1,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_tur,pUnit_calc,pUnit_conv)),1.1*np.max(pressure_conversion(p_boundary_tur,pUnit_calc,pUnit_conv))])\n",
"\n", "\n",
"axs2[1,1].set_title('Velocity Turbine')\n", "axs2[1,1].set_title('Velocity Turbine')\n",
"axs2[1,1].plot(t_vec,v_boundary_tur)\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_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\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", "axs2[1,1].set_ylim([0.95*np.min(v_boundary_tur),1.05*np.max(v_boundary_tur)])\n",
"\n", "\n",
"fig2.tight_layout()\n", "fig2.tight_layout()\n",
"plt.show()" "plt.show()"

288
Regler/Pegelregler_test.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 20, "execution_count": 27,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -23,85 +23,108 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 21, "execution_count": 28,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# define constants\n", "# define constants\n",
"\n", "\n",
"#Turbine\n", " # for physics\n",
"Q_nenn = 0.85 # m³/s\n", "g = 9.81 # [m/s²] gravitational acceleration \n",
"p_nenn = pressure_conversion(10.6,'bar','Pa')\n", "rho = 1000. # [kg/m³] density of water \n",
"closing_time = 480. #s\n", "pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n",
"pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n",
"\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", " # for Turbine\n",
"target_level = 8. # m\n", "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n",
"Kp = 0.01\n", "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n",
"Ti = 200.\n", "Tur_closingTime = 90. # [s] closing time of turbine\n",
"deadband_range = 0.05 # m\n",
"\n", "\n",
"# reservoir\n",
"initial_level = target_level\n",
"initial_influx = Q_nenn/2 # initial influx of volume to the reservoir [m³/s]\n",
"initial_pressure_unit = 'Pa' # DO NOT CHANGE! for pressure conversion in print statements and plot labels \n",
"conversion_pressure_unit = 'bar' # for pressure conversion in print statements and plot labels\n",
"area_base = 74. # total base are of the cuboid reservoir [m²] \n",
"area_outflux = 1. # outflux area of the reservoir, given by pipeline area [m²]\n",
"critical_level_low = 0. # for yet-to-be-implemented warnings[m]\n",
"critical_level_high = np.inf # for yet-to-be-implemented warnings[m]\n",
"\n", "\n",
"p0 = rho*g*initial_level-0.5*rho*(initial_influx/area_outflux)**2\n", " # for PI controller\n",
"Con_targetLevel = 8. # [m]\n",
"Con_K_p = 0.1 # [-] proportional constant of PI controller\n",
"Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\n",
"Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n",
"\n", "\n",
"# offset the pressure in front of the turbine to get realisitc fluxes\n",
"h_fict = 100\n",
"offset_pressure = rho*g*h_fict\n",
"\n", "\n",
"t_max = 1e4 #s\n", " # for pipeline\n",
"dt = 1e-2 # simulation timestep\n", "Pip_length = (535.+478.) # [m] length of pipeline\n",
"nt = int(t_max//dt) # number of simulation steps of reservoir in between timesteps of pipeline \n", "Pip_dia = 0.9 # [m] diameter of pipeline\n",
"Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n",
"Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n",
"Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n",
"Pip_n_seg = 50 # [-] number of pipe segments in discretization\n",
"Pip_f_D = 0.014 # [-] Darcy friction factor\n",
"Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n",
" # derivatives of the pipeline constants\n",
"Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n",
"Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n",
"Pip_nn = Pip_n_seg+1 # [1] number of nodes\n",
"Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n",
"Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n",
"\n", "\n",
"t_vec = np.arange(0,nt+1,1)*dt\n", "\n",
"\n" " # for reservoir\n",
"Res_area_base = 10. # [m²] total base are of the cuboid reservoir \n",
"Res_area_out = Pip_area # [m²] outflux area of the reservoir, given by pipeline area\n",
"Res_level_crit_lo = 0. # [m] for yet-to-be-implemented warnings\n",
"Res_level_crit_hi = np.inf # [m] for yet-to-be-implemented warnings\n",
"Res_dt_approx = 1e-3 # [s] approx. timestep of reservoir time evolution to ensure numerical stability (see Res_nt why approx.)\n",
"Res_nt = max(1,int(Pip_dt//Res_dt_approx)) # [1] number of timesteps of the reservoir time evolution within one timestep of the pipeline\n",
"Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n",
"\n",
" # for general simulation\n",
"flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n",
"level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n",
"simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n",
"nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n",
"t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 22, "execution_count": 29,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# create objects\n", "# create objects\n",
"offset_pressure = pressure_conversion(Pip_head,'mws',pUnit_calc)\n",
"\n", "\n",
"V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,dt)\n", "# Upstream reservoir\n",
"V.set_steady_state(initial_influx,initial_level,conversion_pressure_unit)\n", "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n",
"reservoir.set_steady_state(flux_init,level_init)\n",
"\n", "\n",
"T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time,dt)\n", "# downstream turbine\n",
"T1.set_steady_state(initial_influx,p0+offset_pressure)\n", "turbine = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime,Pip_dt,pUnit_conv)\n",
"turbine.set_steady_state(flux_init,reservoir.get_current_pressure()+offset_pressure)\n",
"\n", "\n",
"Pegelregler = PI_controller_class(target_level,deadband_range,Kp,Ti,dt)" "\n",
"# level controll\n",
"level_control = PI_controller_class(Con_targetLevel,Con_deadbandRange,Con_K_p,Con_T_i,Pip_dt)\n",
"level_control.set_control_variable(turbine.get_current_LA(),display_warning=False)\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 23, "execution_count": 30,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"level_vec = np.full(nt+1,V.level)\n", "level_vec = np.zeros_like(t_vec)\n",
"LA_ist_vec = np.full(nt+1,T1.LA)\n", "level_vec[0] = level_init\n",
"LA_soll_vec = np.full(nt+1,T1.LA)\n", "LA_ist_vec = np.zeros_like(t_vec)\n",
"Q_vec = np.full(nt+1,initial_influx)\n", "LA_ist_vec[0] = turbine.get_current_LA()\n",
"\n", "LA_soll_vec = np.zeros_like(t_vec)\n",
"Pegelregler.control_variable = T1.get_current_LA()" "LA_soll_vec[0] = turbine.get_current_LA()\n",
"Q_vec = np.zeros_like(t_vec)\n",
"Q_vec[0] = turbine.get_current_Q()"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 24, "execution_count": 31,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -109,105 +132,20 @@
"output_type": "stream", "output_type": "stream",
"text": [ "text": [
"0.0\n", "0.0\n",
"100.0\n", "40.52\n",
"200.0\n", "81.04\n",
"300.0\n", "121.56\n",
"400.0\n", "162.08\n",
"500.0\n", "202.6\n",
"600.0\n", "243.12\n",
"700.0\n", "283.64\n",
"800.0\n", "324.16\n",
"900.0\n", "364.68\n",
"1000.0\n", "405.2\n",
"1100.0\n", "445.72\n",
"1200.0\n", "486.24\n",
"1300.0\n", "526.76\n",
"1400.0\n", "567.28\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"
] ]
} }
], ],
@@ -216,39 +154,34 @@
"\n", "\n",
"for i in range(nt+1):\n", "for i in range(nt+1):\n",
"\n", "\n",
" if np.mod(i,1e4) == 0:\n", " if np.mod(i,1e3) == 0:\n",
" print(t_vec[i])\n", " print(t_vec[i])\n",
"\n", "\n",
" if i == 0.2*(nt+1):\n", " if i > 0.1*(nt+1):\n",
" V.set_influx(0.)\n", " reservoir.set_influx(0.)\n",
" elif i == 0.5*(nt+1):\n",
" V.set_influx(0.5*Q_nenn)\n",
" elif i == 0.8*(nt+1):\n",
" V.set_influx(Q_nenn)\n",
"\n", "\n",
"\n", " p = reservoir.get_current_pressure()\n",
" p = V.get_current_pressure()\n", " level_control.update_control_variable(reservoir.level)\n",
" Pegelregler.update_control_variable(V.level)\n", " LA_soll = level_control.get_current_control_variable()\n",
" LA_soll = Pegelregler.get_current_control_variable()\n", " turbine.update_LA(LA_soll)\n",
" T1.update_LA(LA_soll)\n", " turbine.set_pressure(p+offset_pressure)\n",
" T1.set_pressure(p+offset_pressure)\n",
" LA_soll_vec[i] = LA_soll\n", " LA_soll_vec[i] = LA_soll\n",
" LA_ist_vec[i] = T1.get_current_LA()\n", " LA_ist_vec[i] = turbine.get_current_LA()\n",
" Q_vec[i] = T1.get_current_Q()\n", " Q_vec[i] = turbine.get_current_Q()\n",
"\n", "\n",
" \n", " \n",
" V.set_outflux(Q_vec[i])\n", " reservoir.set_outflux(Q_vec[i],display_warning=False)\n",
"\n", "\n",
" V.timestep_reservoir_evolution() \n", " for it_res in range(Res_nt):\n",
" \n", " reservoir.timestep_reservoir_evolution() \n",
" level_vec[i] = V.get_current_level()\n", " level_vec[i] = reservoir.get_current_level()\n",
" \n", " \n",
" " " "
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 25, "execution_count": 32,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -261,12 +194,12 @@
"axs1[0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs1[0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs1[0].set_ylabel(r'$h$ [$\\mathrm{m}$]')\n", "axs1[0].set_ylabel(r'$h$ [$\\mathrm{m}$]')\n",
"axs1[0].plot(t_vec,level_vec)\n", "axs1[0].plot(t_vec,level_vec)\n",
"axs1[0].set_ylim([0*initial_level,1.5*initial_level])\n", "axs1[0].set_ylim([0*level_init,1.5*level_init])\n",
"axs1[1].set_title('Flux')\n", "axs1[1].set_title('Flux')\n",
"axs1[1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs1[1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs1[1].set_ylabel(r'$Q$ [$\\mathrm{m} / \\mathrm{s}^3$]')\n", "axs1[1].set_ylabel(r'$Q$ [$\\mathrm{m} / \\mathrm{s}^3$]')\n",
"axs1[1].plot(t_vec,Q_vec)\n", "axs1[1].plot(t_vec,Q_vec)\n",
"axs1[1].set_ylim([0,2*initial_influx])\n", "axs1[1].set_ylim([0,2*flux_init])\n",
"axs1[2].set_title('LA')\n", "axs1[2].set_title('LA')\n",
"axs1[2].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs1[2].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs1[2].set_ylabel(r'$LA$ [%]')\n", "axs1[2].set_ylabel(r'$LA$ [%]')\n",
@@ -276,27 +209,6 @@
"fig1.tight_layout()\n", "fig1.tight_layout()\n",
"fig1.show()\n" "fig1.show()\n"
] ]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x151d938bc40>]"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig2 = plt.figure()\n",
"plt.plot(t_vec,Pegelregler.get_error_history())"
]
} }
], ],
"metadata": { "metadata": {

6
Regler/Regler_class_file.py
View File

@@ -84,17 +84,17 @@ class PI_controller_class:
# use a list to be able to append more easily - will get converted to np.array when needed # use a list to be able to append more easily - will get converted to np.array when needed
self.error_history = [0] self.error_history = [0]
self.control_variable = -99
self.cv_lower_limit = lower_limit # limits for the controll variable self.cv_lower_limit = lower_limit # limits for the controll variable
self.cv_upper_limit = upper_limit # limits for the controll variable self.cv_upper_limit = upper_limit # limits for the controll variable
# setter # setter
def set_setpoint(self,setpoint): def set_setpoint(self,setpoint):
self.SP = setpoint self.SP = setpoint
def set_control_variable(self,control_variable, display_warning=True): def set_control_variable(self,control_variable, display_warning=True):
if display_warning == True and self.control_variable != -99: if display_warning == True:
print('WARNING! You are setting the control variable of the PI controller manually \ print('WARNING! You are setting the control variable of the PI controller manually \
and are not using the .update_controll_variable() method') and are not using the .update_controll_variable() method')
self.control_variable = control_variable self.control_variable = control_variable

127
Turbinen/Turbinen_class_file.py
View File

@@ -1,8 +1,10 @@
from time import time
import numpy as np import numpy as np
#importing pressure conversion function #importing pressure conversion function
import sys import sys
import os import os
from pyparsing import alphanums
current = os.path.dirname(os.path.realpath(__file__)) current = os.path.dirname(os.path.realpath(__file__))
parent = os.path.dirname(current) parent = os.path.dirname(current)
sys.path.append(parent) sys.path.append(parent)
@@ -11,7 +13,7 @@ from functions.pressure_conversion import pressure_conversion
class Francis_Turbine: class Francis_Turbine:
# units # units
# make sure that units and print units are the same # 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() # units are used to label graphs and disp units are used to have a bearable format when using pythons print()
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}$'
LA_unit = '%' LA_unit = '%'
@@ -20,30 +22,28 @@ class Francis_Turbine:
velocity_unit = r'$\mathrm{m}/\mathrm{s}$' velocity_unit = r'$\mathrm{m}/\mathrm{s}$'
volume_unit = r'$\mathrm{m}^3$' volume_unit = r'$\mathrm{m}^3$'
density_unit_print = 'kg/m³' density_unit_disp = 'kg/m³'
flux_unit_print = 'm³/s' flux_unit_disp = 'm³/s'
LA_unit_print = '%' LA_unit_disp = '%'
pressure_unit_print = 'mWS' time_unit_disp = 's'
time_unit_print = 's' velocity_unit_disp = 'm/s'
velocity_unit_print = 'm/s' volume_unit_disp = 'm³'
volume_unit_print = 'm³'
g = 9.81 # m/s² gravitational acceleration g = 9.81 # m/s² gravitational acceleration
# init # init
def __init__(self, Q_nenn,p_nenn,t_closing=-1.,timestep=-1.): def __init__(self, Q_nenn,p_nenn,t_closing,timestep,pressure_unit_disp):
self.Q_n = Q_nenn # nominal flux self.Q_n = Q_nenn # nominal flux
self.p_n = p_nenn # nominal pressure self.p_n = p_nenn # nominal pressure
self.LA_n = 1. # 100% # nominal Leitapparatöffnung 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.dt = timestep # simulation timestep
self.t_c = t_closing # closing time self.t_c = t_closing # closing time
self.d_LA_max_dt = 1/t_closing # maximal change of LA per second self.d_LA_max_dt = 1/t_closing # maximal change of LA per second
self.pressure_unit_disp = pressure_unit_disp
# initialize for get_info() - parameters will be converted to display -1 if not overwritten # 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.p = pressure_conversion(-1,self.pressure_unit_disp,self.pressure_unit)
self.Q = -1. self.Q = -1.
self.LA = -0.01 self.LA = -0.01
@@ -54,19 +54,22 @@ class Francis_Turbine:
self.LA = LA self.LA = LA
# warn user, that the .set_LA() method should not be used ot set LA manually # warn user, that the .set_LA() method should not be used ot set LA manually
if display_warning == True: if display_warning == True:
print('Consider using the .update_LA() method instead of setting LA manually') print('You are setting the guide vane opening of the turbine manually. \n \
This is not an intended use of this method. \n \
def set_timestep(self,timestep,display_warning=True): Refer to the .update_LA() method instead.')
# 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): def set_pressure(self,pressure):
# set pressure in front of the turbine # set pressure in front of the turbine
self.p = pressure self.p = pressure
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)
if ss_LA < 0 or ss_LA > 1:
raise Exception('LA out of range [0;1]')
self.set_LA(ss_LA,display_warning=False)
#getter #getter
def get_current_Q(self): def get_current_Q(self):
# return the flux through the turbine, based on the current pressure in front # return the flux through the turbine, based on the current pressure in front
@@ -80,10 +83,13 @@ class Francis_Turbine:
def get_current_LA(self): def get_current_LA(self):
return self.LA return self.LA
def get_current_pressure(self):
return pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_disp)
def get_info(self, full = False): def get_info(self, full = False):
new_line = '\n' new_line = '\n'
p = pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_print) p = pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_disp)
p_n = pressure_conversion(self.p_n,self.pressure_unit,self.pressure_unit_print) p_n = pressure_conversion(self.p_n,self.pressure_unit,self.pressure_unit_disp)
if full == True: if full == True:
@@ -91,33 +97,34 @@ class Francis_Turbine:
print_str = (f"Turbine has the following attributes: {new_line}" print_str = (f"Turbine has the following attributes: {new_line}"
f"----------------------------- {new_line}" f"----------------------------- {new_line}"
f"Type = Francis {new_line}" f"Type = Francis {new_line}"
f"Nominal flux = {self.Q_n:<10} {self.flux_unit_print} {new_line}" f"Nominal flux = {self.Q_n:<10} {self.flux_unit_disp} {new_line}"
f"Nominal pressure = {round(p_n,3):<10} {self.pressure_unit_print}{new_line}" f"Nominal pressure = {round(p_n,3):<10} {self.pressure_unit_disp}{new_line}"
f"Nominal LA = {self.LA_n*100:<10} {self.LA_unit_print} {new_line}" f"Nominal LA = {self.LA_n*100:<10} {self.LA_unit_disp} {new_line}"
f"Closing time = {self.t_c:<10} {self.time_unit_print} {new_line}" f"Closing time = {self.t_c:<10} {self.time_unit_disp} {new_line}"
f"Current flux = {self.Q:<10} {self.flux_unit_print} {new_line}" f"Current flux = {self.Q:<10} {self.flux_unit_disp} {new_line}"
f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}" f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_disp} {new_line}"
f"Current LA = {self.LA*100:<10} {self.LA_unit_print} {new_line}" f"Current LA = {self.LA*100:<10} {self.LA_unit_disp} {new_line}"
f"Simulation timestep = {self.dt:<10} {self.time_unit_print} {new_line}" f"Simulation timestep = {self.dt:<10} {self.time_unit_disp} {new_line}"
f"----------------------------- {new_line}") f"----------------------------- {new_line}")
else: else:
# :<10 pads the self.value to be 10 characters wide # :<10 pads the self.value to be 10 characters wide
print_str = (f"The current attributes are: {new_line}" print_str = (f"The current attributes are: {new_line}"
f"----------------------------- {new_line}" f"----------------------------- {new_line}"
f"Current flux = {self.Q:<10} {self.flux_unit_print} {new_line}" f"Current flux = {self.Q:<10} {self.flux_unit_disp} {new_line}"
f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}" f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_disp} {new_line}"
f"Current LA = {self.LA*100:<10} {self.LA_unit_print} {new_line}" f"Current LA = {self.LA*100:<10} {self.LA_unit_disp} {new_line}"
f"----------------------------- {new_line}") f"----------------------------- {new_line}")
print(print_str) print(print_str)
# methods # update methods
def update_LA(self,LA_soll): def update_LA(self,LA_soll):
# update the Leitappartöffnung and consider the restrictions of the closing time of the turbine # 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 = 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_max = self.d_LA_max_dt*self.dt # calculate the maximum possible change in LA based on the given timestep
LA_diff = np.sign(LA_diff)*np.min(np.abs([LA_diff,LA_diff_max])) # calulate the correct change in LA LA_diff = np.sign(LA_diff)*np.min(np.abs([LA_diff,LA_diff_max])) # calulate the correct change in LA
# make sure that the LA is not out of the range [0;1]
LA_new = self.LA-LA_diff LA_new = self.LA-LA_diff
if LA_new < 0.: if LA_new < 0.:
LA_new = 0. LA_new = 0.
@@ -125,10 +132,42 @@ class Francis_Turbine:
LA_new = 1. LA_new = 1.
self.set_LA(LA_new,display_warning=False) self.set_LA(LA_new,display_warning=False)
def set_steady_state(self,ss_flux,ss_pressure): # methods
# calculate and set steady state LA, that allows the flow of ss_flux at ss_pressure through the def converge(self,convergence_parameters):
# turbine at the steady state LA # small numerical disturbances (~1e-12 m/s) in the velocity can get amplified at the turbine node, because the new velocity of the turbine and the
ss_LA = self.LA_n*ss_flux/self.Q_n*np.sqrt(self.p_n/ss_pressure) # new pressure from the forward characteristic are not compatible.
if ss_LA < 0 or ss_LA > 1: eps = 1e-12 # convergence criterion: iteration change < eps
raise Exception('LA out of range [0;1]') iteration_change = 1. # change in Q from one iteration to the next
self.set_LA(ss_LA,display_warning=False) i = 0 # safety variable. break loop if it exceeds 1e6 iterations
g = self.g # gravitational acceleration
p = convergence_parameters[0] # pressure at second to last node (see method of characterisctics - boundary condidtions)
v = convergence_parameters[1] # velocity at second to last node (see method of characterisctics - boundary condidtions)
D = convergence_parameters[2] # diameter of the pipeline
area_pipe = convergence_parameters[3] # area of the pipeline
alpha = convergence_parameters[4] # elevation angle of the pipeline
f_D = convergence_parameters[5] # Darcy friction coefficient
c = convergence_parameters[6] # pressure wave propagtation velocity
rho = convergence_parameters[7] # density of the liquid
dt = convergence_parameters[8] # timestep of the characteristic method
p_old = self.get_current_pressure()
Q_old = self.get_current_Q()
v_old = Q_old/area_pipe
while iteration_change > eps:
self.set_pressure(p_old)
Q_new = self.get_current_Q()
v_new = Q_new/area_pipe
p_new = p-rho*c*(v_old-v)+rho*c*dt*g*np.sin(alpha)-f_D*rho*c*dt/(2*D)*abs(v)*v
iteration_change = abs(Q_old-Q_new)
Q_old = Q_new.copy()
p_old = p_new.copy()
v_old = v_new.copy()
i = i+1
if i == 1e6:
print('did not converge')
break
# print(i)
self.Q = Q_new

370
Turbinen/Turbinen_test_steady_state.ipynb Normal file
View File

@@ -0,0 +1,370 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from Turbinen_class_file import Francis_Turbine\n",
"\n",
"import sys\n",
"import os\n",
"current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n",
"parent = os.path.dirname(current)\n",
"sys.path.append(parent)\n",
"from functions.pressure_conversion import pressure_conversion\n",
"from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class\n",
"from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class\n",
"from Regler.Regler_class_file import PI_controller_class"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"# define constants\n",
"\n",
" # for physics\n",
"g = 9.81 # [m/s²] gravitational acceleration \n",
"rho = 1000. # [kg/m³] density of water \n",
"pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n",
"pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n",
"\n",
"\n",
" # for Turbine\n",
"Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n",
"Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n",
"Tur_closingTime = 90. # [s] closing time of turbine\n",
"\n",
"\n",
" # for PI controller\n",
"Con_targetLevel = 8. # [m]\n",
"Con_K_p = 0.1 # [-] proportional constant of PI controller\n",
"Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\n",
"Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n",
"\n",
"\n",
" # for pipeline\n",
"Pip_length = (535.+478.) # [m] length of pipeline\n",
"Pip_dia = 0.9 # [m] diameter of pipeline\n",
"Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n",
"Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n",
"Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n",
"Pip_n_seg = 50 # [-] number of pipe segments in discretization\n",
"Pip_f_D = 0.014 # [-] Darcy friction factor\n",
"Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n",
" # derivatives of the pipeline constants\n",
"Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n",
"Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n",
"Pip_nn = Pip_n_seg+1 # [1] number of nodes\n",
"Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n",
"Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n",
"\n",
"\n",
" # for reservoir\n",
"Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n",
"Res_area_out = Pip_area # [m²] outflux area of the reservoir, given by pipeline area\n",
"Res_level_crit_lo = 0. # [m] for yet-to-be-implemented warnings\n",
"Res_level_crit_hi = np.inf # [m] for yet-to-be-implemented warnings\n",
"Res_dt_approx = 1e-3 # [s] approx. timestep of reservoir time evolution to ensure numerical stability (see Res_nt why approx.)\n",
"Res_nt = max(1,int(Pip_dt//Res_dt_approx)) # [1] number of timesteps of the reservoir time evolution within one timestep of the pipeline\n",
"Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n",
"\n",
" # for general simulation\n",
"flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n",
"level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n",
"simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n",
"nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n",
"t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"# create objects\n",
"\n",
"# Upstream reservoir\n",
"reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n",
"reservoir.set_steady_state(flux_init,level_init)\n",
"\n",
"# pipeline\n",
"pipe = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_n_seg,Pip_angle,Pip_f_D,Pip_pw_vel,Pip_dt,pUnit_conv,rho)\n",
"pipe.set_steady_state(flux_init,level_init,Res_area_base,Pip_x_vec,Pip_h_vec)\n",
"\n",
"# downstream turbine\n",
"turbine = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime,Pip_dt,pUnit_conv)\n",
"turbine.set_steady_state(flux_init,pipe.get_current_pressure_distribution()[-1])\n",
"\n",
"# influx setting turbine\n",
"turbine_in = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime/2,Pip_dt,pUnit_conv)\n",
"turbine_in.set_steady_state(flux_init,Tur_p_nenn)\n",
"\n",
"# level controll\n",
"level_control = PI_controller_class(Con_targetLevel,Con_deadbandRange,Con_K_p,Con_T_i,Pip_dt)\n",
"level_control.set_control_variable(turbine.get_current_LA(),display_warning=False)\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# initialization for Timeloop\n",
"\n",
"v_old = pipe.get_current_velocity_distribution()\n",
"v_min = pipe.get_current_velocity_distribution()\n",
"v_max = pipe.get_current_velocity_distribution()\n",
"Q_old = pipe.get_current_flux_distribution()\n",
"Q_min = pipe.get_current_flux_distribution()\n",
"Q_max = pipe.get_current_flux_distribution()\n",
"p_old = pipe.get_current_pressure_distribution()\n",
"p_min = pipe.get_current_pressure_distribution()\n",
"p_max = pipe.get_current_pressure_distribution()\n",
"\n",
"Q_in_vec = np.zeros_like(t_vec)\n",
"Q_in_vec[0] = flux_init\n",
"\n",
"v_boundary_res = np.zeros_like(t_vec)\n",
"v_boundary_tur = np.zeros_like(t_vec)\n",
"Q_boundary_res = np.zeros_like(t_vec)\n",
"Q_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",
"level_vec = np.full_like(t_vec,level_init) # level at the end of each pipeline timestep\n",
"volume_vec = np.full_like(t_vec,reservoir.get_current_volume()) # volume at the end of each pipeline timestep\n",
"\n",
"v_boundary_res[0] = v_old[0]\n",
"v_boundary_tur[0] = v_old[-1] \n",
"Q_boundary_res[0] = Q_old[0]\n",
"Q_boundary_tur[0] = Q_old[-1]\n",
"p_boundary_res[0] = p_old[0]\n",
"p_boundary_tur[0] = p_old[-1]\n",
"\n",
"LA_soll_vec = np.full_like(t_vec,turbine.get_current_LA())\n",
"LA_ist_vec = np.full_like(t_vec,turbine.get_current_LA())\n",
"\n",
"LA_soll_vec2 = np.full_like(t_vec,turbine_in.get_current_LA())\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib qt5\n",
"# Con_T_ime loop\n",
"\n",
"# create a figure and subplots to display the velocity and pressure distribution across the pipeline in each pipeline step\n",
"fig1,axs1 = plt.subplots(2,1)\n",
"fig1.suptitle(str(0) +' s / '+str(round(t_vec[-1],2)) + ' s' )\n",
"axs1[0].set_title('Pressure distribution in pipeline')\n",
"axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs1[0].set_ylabel(r'$p$ ['+pUnit_conv+']')\n",
"axs1[1].set_title('Flux distribution in pipeline')\n",
"axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs1[1].set_ylabel(r'$Q$ [$\\mathrm{m}^3 / \\mathrm{s}$]')\n",
"lo_p, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,pUnit_calc, pUnit_conv),marker='.')\n",
"lo_q, = axs1[1].plot(Pip_x_vec,Q_old,marker='.')\n",
"lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n",
"lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n",
"lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n",
"lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n",
"\n",
"axs1[0].autoscale()\n",
"axs1[1].autoscale()\n",
"\n",
"fig1.tight_layout()\n",
"fig1.show()\n",
"plt.pause(1)\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"convergence_parameters = [p_old[-2],v_old[-2],Pip_dia,Pip_area,Pip_angle,Pip_f_D,Pip_pw_vel,rho,Pip_dt]\n",
"\n",
"# loop through Con_T_ime steps of the pipeline\n",
"for it_pipe in range(1,nt+1):\n",
"\n",
" turbine_in.update_LA(LA_soll_vec2[it_pipe])\n",
" turbine_in.set_pressure(Tur_p_nenn)\n",
" Q_in_vec[it_pipe] = turbine_in.get_current_Q()\n",
" reservoir.set_influx(Q_in_vec[it_pipe])\n",
"\n",
"# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n",
" # set initial condition for the reservoir Con_T_ime evolution calculted with e-RK4\n",
" reservoir.set_pressure(p_old[0],display_warning=False)\n",
" reservoir.set_outflux(Q_old[0],display_warning=False)\n",
" # calculate the Con_T_ime evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n",
" for it_res in range(Res_nt):\n",
" reservoir.timestep_reservoir_evolution() \n",
" level_vec[it_pipe] = reservoir.get_current_level() \n",
" volume_vec[it_pipe] = reservoir.get_current_volume() \n",
"\n",
" # get the control variable\n",
" level_control.update_control_variable(level_vec[it_pipe])\n",
" LA_soll_vec[it_pipe] = level_control.get_current_control_variable()\n",
" \n",
" # change the Leitapparatöffnung based on the target value\n",
" turbine.update_LA(LA_soll_vec[it_pipe])\n",
" LA_ist_vec[it_pipe] = turbine.get_current_LA()\n",
"\n",
" # set boundary condition for the next timestep of the characterisCon_T_ic method\n",
" turbine.set_pressure(p_old[-1])\n",
" convergence_parameters[0] = p_old[-2]\n",
" convergence_parameters[1] = v_old[-2]\n",
" turbine.converge(convergence_parameters)\n",
" p_boundary_res[it_pipe] = reservoir.get_current_pressure()\n",
" v_boundary_tur[it_pipe] = 1/Pip_area*turbine.get_current_Q()\n",
" Q_boundary_tur[it_pipe] = turbine.get_current_Q()\n",
"\n",
" # the the boundary condition 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",
" pipe.v[0] = (0.8*pipe.v[0]+0.2*reservoir.get_current_outflux()/Res_area_out)\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",
" Q_boundary_res[it_pipe] = pipe.get_current_flux_distribution()[0]\n",
"\n",
" # perform the next timestep via the characterisCon_T_ic 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",
" Q_old = pipe.get_current_flux_distribution()\n",
"\n",
"\n",
" # plot some stuff\n",
" # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n",
" lo_p.remove()\n",
" lo_pmin.remove()\n",
" lo_pmax.remove()\n",
" lo_q.remove()\n",
" lo_qmin.remove()\n",
" lo_qmax.remove()\n",
" # plot new pressure and velocity distribution in the pipeline\n",
" lo_p, = axs1[0].plot(Pip_x_vec,pipe.get_current_pressure_distribution(disp=True),marker='.',c='blue')\n",
" lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n",
" lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n",
" lo_q, = axs1[1].plot(Pip_x_vec,pipe.get_current_flux_distribution(),marker='.',c='blue')\n",
" lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n",
" lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\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.001) "
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# plot Con_T_ime evolution of boundary pressure and velocity as well as the reservoir level\n",
"\n",
"fig2,axs2 = plt.subplots(1,1)\n",
"axs2.set_title('Level and Volume reservoir')\n",
"axs2.plot(t_vec,level_vec,label='level')\n",
"axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2.set_ylabel(r'$h$ [m]')\n",
"x_twin_00 = axs2.twinx()\n",
"x_twin_00.set_ylabel(r'$V$ [$\\mathrm{m}^3$]')\n",
"x_twin_00.plot(t_vec,volume_vec)\n",
"axs2.legend()\n",
"\n",
"fig2,axs2 = plt.subplots(1,1)\n",
"axs2.set_title('LA')\n",
"axs2.plot(t_vec,100*LA_soll_vec,label='Target')\n",
"axs2.plot(t_vec,100*LA_ist_vec,label='Actual')\n",
"axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2.set_ylabel(r'$LA$ [%]')\n",
"axs2.legend()\n",
"\n",
"fig2,axs2 = plt.subplots(1,1)\n",
"axs2.set_title('Pressure reservoir and turbine')\n",
"axs2.plot(t_vec,pressure_conversion(p_boundary_res,pUnit_calc, pUnit_conv),label='Reservoir')\n",
"axs2.plot(t_vec,pressure_conversion(p_boundary_tur,pUnit_calc, pUnit_conv),label='Turbine')\n",
"axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n",
"axs2.legend()\n",
"\n",
"fig2,axs2 = plt.subplots(1,1)\n",
"axs2.set_title('Fluxes')\n",
"axs2.plot(t_vec,Q_boundary_res,label='Outflux')\n",
"axs2.plot(t_vec,Q_in_vec,label='Influx')\n",
"axs2.plot(t_vec,Q_boundary_tur,label='Flux Turbine')\n",
"axs2.set_ylim(-2*Tur_Q_nenn,+2*Tur_Q_nenn)\n",
"axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2.set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n",
"axs2.legend()\n",
"\n",
"fig2,axs2 = plt.subplots(1,1)\n",
"axs2.set_title('Min and Max Pressure')\n",
"axs2.plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n",
"axs2.plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n",
"axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n",
"\n",
"fig2,axs2 = plt.subplots(1,1)\n",
"axs2.set_title('Min and Max Fluxes')\n",
"axs2.plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n",
"axs2.plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n",
"axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs2.set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n",
"\n",
"# axs2[0,1].legend()\n",
"# axs2[1,0].legend()\n",
"# axs2[1,1].legend()\n",
"# # axs2[2,0].legend()\n",
"# # axs2[2,1].legend()\n",
"\n",
"\n",
"fig2.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
}

168
Turbinen/old/convergence_turbine.py Normal file
View File

@@ -0,0 +1,168 @@
import numpy as np
#importing pressure conversion function
import sys
import os
current = os.path.dirname(os.path.realpath(__file__))
parent = os.path.dirname(current)
sys.path.append(parent)
from functions.pressure_conversion import pressure_conversion
class Francis_Turbine_test:
# 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$'
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 = timestep
# warn user, that the .set_LA() method should not be used ot set LA manually
if display_warning == True:
print('WARNING: You are changing the timestep of the turbine simulation. This has implications on the simulated closing speed!')
def set_pressure(self,pressure):
# set pressure in front of the turbine
self.p = pressure
#getter
def get_current_Q(self):
# return the flux through the turbine, based on the current pressure in front
# of the turbine and the Leitapparatöffnung
if self.p < 0:
self.Q = 0
else:
self.Q = self.Q_n*(self.LA/self.LA_n)*np.sqrt(self.p/self.p_n)
return self.Q
def get_current_pressure(self):
return self.p
def get_current_LA(self):
return self.LA
def get_info(self, full = False):
new_line = '\n'
p = pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_print)
p_n = pressure_conversion(self.p_n,self.pressure_unit,self.pressure_unit_print)
if full == True:
# :<10 pads the self.value to be 10 characters wide
print_str = (f"Turbine has the following attributes: {new_line}"
f"----------------------------- {new_line}"
f"Type = Francis {new_line}"
f"Nominal flux = {self.Q_n:<10} {self.flux_unit_print} {new_line}"
f"Nominal pressure = {round(p_n,3):<10} {self.pressure_unit_print}{new_line}"
f"Nominal LA = {self.LA_n*100:<10} {self.LA_unit_print} {new_line}"
f"Closing time = {self.t_c:<10} {self.time_unit_print} {new_line}"
f"Current flux = {self.Q:<10} {self.flux_unit_print} {new_line}"
f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}"
f"Current LA = {self.LA*100:<10} {self.LA_unit_print} {new_line}"
f"Simulation timestep = {self.dt:<10} {self.time_unit_print} {new_line}"
f"----------------------------- {new_line}")
else:
# :<10 pads the self.value to be 10 characters wide
print_str = (f"The current attributes are: {new_line}"
f"----------------------------- {new_line}"
f"Current flux = {self.Q:<10} {self.flux_unit_print} {new_line}"
f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}"
f"Current LA = {self.LA*100:<10} {self.LA_unit_print} {new_line}"
f"----------------------------- {new_line}")
print(print_str)
# methods
def update_LA(self,LA_soll):
# update the Leitappartöffnung and consider the restrictions of the closing time of the turbine
LA_diff = self.LA-LA_soll # calculate the difference to the target LA
LA_diff_max = self.d_LA_max_dt*self.dt # calculate the maximum change in LA based on the given timestep
LA_diff = np.sign(LA_diff)*np.min(np.abs([LA_diff,LA_diff_max])) # calulate the correct change in LA
LA_new = self.LA-LA_diff
if LA_new < 0.:
LA_new = 0.
elif LA_new > 1.:
LA_new = 1.
self.set_LA(LA_new,display_warning=False)
def set_steady_state(self,ss_flux,ss_pressure):
# calculate and set steady state LA, that allows the flow of ss_flux at ss_pressure through the
# turbine at the steady state LA
ss_LA = self.LA_n*ss_flux/self.Q_n*np.sqrt(self.p_n/ss_pressure)
if ss_LA < 0 or ss_LA > 1:
raise Exception('LA out of range [0;1]')
self.set_LA(ss_LA,display_warning=False)
def converge(self,area_pipe,pressure_s2l_node,velocity_s2l_node,alpha,f_D,dt):
eps = 1e-9
error = 1.
i = 0
p = pressure_s2l_node
v = velocity_s2l_node
rho = 1000
g = self.g
c = 400
D = area_pipe
p_old = self.get_current_pressure()
Q_old = self.get_current_Q()
v_old = Q_old/area_pipe
while error > eps:
self.set_pressure(p_old)
Q_new = self.get_current_Q()
v_new = Q_new/area_pipe
p_new = p-rho*c*(v_old-v)+rho*c*dt*g*np.sin(alpha)-f_D*rho*c*dt/(2*D)*abs(v)*v
error = abs(Q_old-Q_new)
Q_old = Q_new.copy()
p_old = p_new.copy()
v_old = v_new.copy()
i = i+1
if i == 1e6:
print('did not converge')
break
self.Q = Q_new

222
Untertweng.ipynb → Turbinen/old/turbine_convergence_test.ipynb
View File

@@ -8,11 +8,18 @@
"source": [ "source": [
"import numpy as np\n", "import numpy as np\n",
"import matplotlib.pyplot as plt\n", "import matplotlib.pyplot as plt\n",
"from convergence_turbine import Francis_Turbine_test\n",
"\n", "\n",
"#importing pressure conversion function\n",
"import sys\n",
"import os\n",
"current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n",
"parent = os.path.dirname(current)\n",
"sys.path.append(parent)\n",
"from functions.pressure_conversion import pressure_conversion\n", "from functions.pressure_conversion import pressure_conversion\n",
"from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class\n", "from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class\n",
"from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class\n", "from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class\n",
"from Turbinen.Turbinen_class_file import Francis_Turbine" "from Regler.Regler_class_file import PI_controller_class"
] ]
}, },
{ {
@@ -21,55 +28,59 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"#define constants\n", "%matplotlib qt5\n",
"\n", "\n",
"#Turbine\n", "#Turbine\n",
"Q_nenn = 0.85 # m³/s\n", "Q_nenn = 0.85 # m³/s\n",
"p_nenn = pressure_conversion(10.6,'bar','Pa')\n", "p_nenn = pressure_conversion(10.6,'bar','Pa')\n",
"closing_time = 5 #s\n", "closing_time = 90. #s\n",
"\n",
"#define constants pipe\n",
"\n", "\n",
"# physics\n",
"g = 9.81 # gravitational acceleration [m/s²]\n", "g = 9.81 # gravitational acceleration [m/s²]\n",
"rho = 1000. # density of water [kg/m³]\n", "rho = 1000. # density of water [kg/m³]\n",
"\n", "\n",
"# pipeline\n", "# define controller constants\n",
"target_level = 8. # m\n",
"Kp = 0.1\n",
"Ti = 1000.\n",
"deadband_range = 0.05 # m\n",
"\n",
"L = 535.+478. # length of pipeline [m]\n", "L = 535.+478. # length of pipeline [m]\n",
"D = 0.9 # pipe diameter [m]\n", "D = 0.9 # pipe diameter [m]\n",
"A_pipe = D**2/4*np.pi # pipeline area\n", "h_res = target_level # water level in upstream reservoir [m]\n",
"h_pipe = 105 # hydraulic head without reservoir [m] \n", "n = 50 # number of pipe segments in discretization\n",
"alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n", "nt = 10000 # number of time steps after initial conditions\n",
"n = 50 # number of pipe segments in discretization # initial flow velocity [m/s]\n",
"f_D = 0.014 # Darcy friction factor\n", "f_D = 0.014 # Darcy friction factor\n",
"c = 400. # propagation velocity of the pressure wave [m/s]\n", "c = 400. # propagation velocity of the pressure wave [m/s]\n",
"# consider prescribing a total simulation time and deducting the number of timesteps from that\n", "h_pipe = 105. # hydraulic head without reservoir [m] \n",
"nt = 9000 # number of time steps after initial conditions\n", "alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n",
"\n", "\n",
"# derivatives of the pipeline constants\n", "# define constants reservoir\n",
"conversion_pressure_unit = 'mWS'\n",
"\n",
"# preparing the discretization and initial conditions\n",
"initial_flux = Q_nenn/1.1 # m³/s\n",
"initial_level = h_res # m\n",
"dx = L/n # length of each pipe segment\n", "dx = L/n # length of each pipe segment\n",
"dt = dx/c # timestep according to method of characterisitics\n", "dt = dx/c # timestep according to method of characterisitics\n",
"nn = n+1 # number of nodes\n", "nn = n+1 # number of nodes\n",
"initial_level = 8. # 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", "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", "t_vec = np.arange(0,nt,1)*dt # time vector\n",
"h_vec = np.arange(0,nn,1)*h_pipe/n # hydraulic head of pipeline at each node\n", "h_vec = np.arange(0,nn,1)*h_pipe/n # hydraulic head of pipeline at each node\n",
"\n", "\n",
"\n", "\n",
"# define constants reservoir\n",
"conversion_pressure_unit = 'mWS'\n",
"\n", "\n",
"# reservoir\n", "area_base = 75. # m²\n",
"# replace influx by vector\n", "area_pipe = (D/2)**2*np.pi # m²\n",
"initial_flux = Q_nenn/1.1 # initial influx of volume to the reservoir [m³/s]\n", "critical_level_low = 0. # m\n",
"initial_pressure_unit = 'Pa' # DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", "critical_level_high = 100. # m\n",
"conversion_pressure_unit = 'bar' # for pressure conversion in print statements and plot labels\n",
"area_base = 74. # total base are of the cuboid reservoir [m²] \n",
"area_outflux = A_pipe # outlfux area of the reservoir, given by pipeline area [m²]\n",
"critical_level_low = 0. # for yet-to-be-implemented warnings[m]\n",
"critical_level_high = np.inf # for yet-to-be-implemented warnings[m]\n",
"\n", "\n",
"# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n", "# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n",
"nt_eRK4 = 100 # number of simulation steps of reservoir in between timesteps of pipeline \n", "nt_eRK4 = 100 # number of simulation steps of reservoir in between timesteps of pipeline \n",
"simulation_timestep = dt/nt_eRK4\n", "simulation_timestep = dt/nt_eRK4"
"\n",
"\n"
] ]
}, },
{ {
@@ -78,21 +89,24 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# create objects\n", "V = Ausgleichsbecken_class(area_base,area_pipe,critical_level_low,critical_level_high,simulation_timestep)\n",
"\n",
"V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,simulation_timestep)\n",
"V.set_steady_state(initial_flux,initial_level,conversion_pressure_unit)\n", "V.set_steady_state(initial_flux,initial_level,conversion_pressure_unit)\n",
"\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_flux,initial_level,area_base,pl_vec,h_vec)\n", "pipe.set_steady_state(initial_flux,initial_level,area_base,pl_vec,h_vec)\n",
"\n", "\n",
"initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n",
"\n", "\n",
"T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time,timestep=dt)\n", "initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n",
"T1.set_steady_state(initial_flux,initial_pressure_turbine)\n" "T1 = Francis_Turbine_test(Q_nenn,p_nenn,closing_time,timestep=dt)\n",
"T1.set_steady_state(initial_flux,initial_pressure_turbine)\n",
"\n",
"T_in = Francis_Turbine_test(Q_nenn,p_nenn,closing_time/2,timestep=dt)\n",
"T_in.set_steady_state(initial_flux,p_nenn)\n",
"\n",
"Pegelregler = PI_controller_class(target_level,deadband_range,Kp,Ti,dt)\n",
"Pegelregler.control_variable = T1.get_current_LA()"
] ]
}, },
{ {
@@ -103,33 +117,36 @@
"source": [ "source": [
"# initialization for timeloop\n", "# initialization for timeloop\n",
"\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", "# 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", "v_old = pipe.get_current_velocity_distribution()\n",
"p_old = pipe.get_current_pressure_distribution()\n", "p_old = pipe.get_current_pressure_distribution()\n",
"\n", "\n",
"# prepare the vectors in which the temporal evolution of the boundary conditions are stored\n", "# prepare the vectors in which the temporal evolution of the boundary conditions are stored\n",
" # keep in mind, that the velocity at the turbine and the pressure at the reservoir follow from boundary conditions\n", " # keep in mind, that the velocity at the turbine and the pressure at the reservoir are set manually and\n",
" # reservoir level and flow through turbine\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.zeros_like(t_vec)\n", "v_boundary_res = np.zeros_like(t_vec)\n",
"v_boundary_tur = 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_res = np.zeros_like(t_vec)\n",
"p_boundary_tur = np.zeros_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",
"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", "# set the boundary conditions for the first timestep\n",
"v_boundary_res[0] = v_old[0]\n", "v_boundary_res[0] = v_old[0]\n",
"v_boundary_tur[0] = v_old[-1] \n", "v_boundary_tur[0] = v_old[-1] \n",
"p_boundary_res[0] = p_old[0]\n", "p_boundary_res[0] = p_old[0]\n",
"p_boundary_tur[0] = p_old[-1]\n", "p_boundary_tur[0] = p_old[-1]\n",
"\n", "\n",
"LA_soll_vec = np.full_like(t_vec,T1.LA)\n", "LA_soll_vec = np.full_like(t_vec,T1.get_current_LA())\n",
"# LA_soll_vec[500:]= 0\n", "LA_ist_vec = np.full_like(t_vec,T1.get_current_LA())\n",
"LA_ist_vec = np.full_like(t_vec,T1.LA)\n",
"\n", "\n",
"\n" "LA_soll_vec2 = np.full_like(t_vec,T_in.get_current_LA())\n",
"LA_soll_vec2[500:1000] = 0.\n",
"LA_soll_vec2[1000:1500] = 1. \n",
"LA_soll_vec2[1500:2000] = 0.\n",
"LA_soll_vec2[2000:2500] = 0.5 "
] ]
}, },
{ {
@@ -138,26 +155,22 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"%matplotlib qt5\n",
"# time loop\n",
"\n",
"# create a figure and subplots to display the velocity and pressure distribution across the pipeline in each pipeline step\n",
"fig1,axs1 = plt.subplots(2,1)\n", "fig1,axs1 = plt.subplots(2,1)\n",
"fig1.suptitle(str(0) +' s / '+str(round(t_vec[-1],2)) + ' s' )\n",
"axs1[0].set_title('Pressure distribution in pipeline')\n", "axs1[0].set_title('Pressure distribution in pipeline')\n",
"axs1[1].set_title('Velocity distribution in pipeline')\n",
"axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs1[0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\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",
"\n",
"axs1[1].set_title('Velocity distribution in pipeline')\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$ [m/s]')\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", "lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n",
"\n",
"axs1[0].autoscale()\n", "axs1[0].autoscale()\n",
"axs1[1].autoscale()\n", "axs1[1].autoscale()\n",
"\n", "\n",
"fig1.tight_layout()\n", "fig1.tight_layout()\n",
"fig1.show()\n", "plt.pause(1)"
"plt.pause(1)\n"
] ]
}, },
{ {
@@ -166,33 +179,39 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# loop through time steps of the pipeline\n",
"for it_pipe in range(1,pipe.nt+1):\n",
"\n",
" # if it_pipe == 250:\n",
" # V.set_influx(0.)\n",
"\n", "\n",
"for it_pipe in range(1,nt):\n",
"# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n", "# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n",
" \n",
" T_in.update_LA(LA_soll_vec2[it_pipe])\n",
" T_in.set_pressure(p_nenn)\n",
" V.set_influx(T_in.get_current_Q())\n",
"\n",
" # set initial conditions for the reservoir time evolution calculted with e-RK4\n", " # set initial conditions for the reservoir time evolution calculted with e-RK4\n",
" V.set_pressure(p_old[0])\n", " V.set_pressure(p_old[0])\n",
" V.set_outflux(v_old[0]*area_outflux)\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", " # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n",
" for it_res in range(nt_eRK4):\n", " for it_res in range(nt_eRK4):\n",
" V.timestep_reservoir_evolution() \n", " V.timestep_reservoir_evolution() \n",
" level_vec[it_pipe] = V.get_current_level() \n", " level_vec[it_pipe] = V.get_current_level() \n",
"\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", " # change the Leitapparatöffnung based on the target value\n",
" T1.update_LA(LA_soll_vec[it_pipe])\n", " T1.update_LA(LA_soll_vec[it_pipe])\n",
" T1.set_pressure(p_old[-1])\n",
"\n",
" LA_ist_vec[it_pipe] = T1.get_current_LA()\n", " LA_ist_vec[it_pipe] = T1.get_current_LA()\n",
"\n", "\n",
" # set boundary conditions for the next timestep of the characteristic method\n", " # set boundary conditions for the next timestep of the characteristic method\n",
" T1.set_pressure(p_old[-1])\n",
" T1.converge(area_pipe,p_old[-2],v_old[-2],alpha,f_D,dt)\n",
" p_boundary_res[it_pipe] = V.get_current_pressure()\n", " p_boundary_res[it_pipe] = V.get_current_pressure()\n",
" v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_current_Q()\n", " v_boundary_tur[it_pipe] = T1.get_current_Q()/area_pipe\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(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",
" pipe.v[0] = (0.8*pipe.v[0]+0.2*V.get_current_outflux()/area_pipe)\n",
" p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\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", " v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n",
"\n", "\n",
@@ -209,34 +228,13 @@
" 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(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.',c='blue')\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", " lo_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n",
" # lo_02, = axs1[2].plot(level_vec_2,c='blue')\n", " \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)) + '/' + str(round(t_vec[-1],2)))\n",
" fig1.canvas.draw()\n", " fig1.canvas.draw()\n",
" fig1.tight_layout()\n", " fig1.tight_layout()\n",
" fig1.show()\n", " plt.pause(0.000001)"
" plt.pause(0.001) \n",
"\n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.7727272727272726\n"
]
}
],
"source": [
"print(V.get_current_influx())"
] ]
}, },
{ {
@@ -249,7 +247,7 @@
"\n", "\n",
"fig2,axs2 = plt.subplots(3,2)\n", "fig2,axs2 = plt.subplots(3,2)\n",
"axs2[0,0].set_title('Pressure reservoir')\n", "axs2[0,0].set_title('Pressure reservoir')\n",
"axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,initial_pressure_unit, conversion_pressure_unit))\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_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[0,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", "axs2[0,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
"\n", "\n",
@@ -260,7 +258,7 @@
"axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", "axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
"\n", "\n",
"axs2[1,0].set_title('Pressure turbine')\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].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_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[1,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", "axs2[1,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
"\n", "\n",
@@ -282,52 +280,6 @@
"fig2.tight_layout()\n", "fig2.tight_layout()\n",
"plt.show()" "plt.show()"
] ]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The cuboid reservoir has the following attributes: \n",
"----------------------------- \n",
"Base area = 74.0 m² \n",
"Outflux area = 0.636 m² \n",
"Current level = 8.020201323491174 m\n",
"Critical level low = 0.0 m \n",
"Critical level high = inf m \n",
"Volume in reservoir = -- m³ \n",
"Current influx = 0.7727272727272726 m³/s \n",
"Current outflux = -3.529572535198196 m³/s \n",
"Current outflux vel = -5.548 m/s \n",
"Current pipe pressure = 0.479 bar \n",
"Simulation timestep = 0.0005065 s \n",
"Density of liquid = 1000 kg/m³ \n",
"----------------------------- \n",
"\n",
"9.47597527926513\n",
"10.931749235039087\n",
"12.387523190813043\n",
"13.843297146587\n",
"15.299071102360957\n",
"16.75484505813491\n",
"18.210619013908868\n",
"19.666392969682825\n",
"21.12216692545678\n",
"22.577940881230738\n"
]
}
],
"source": [
"V.get_info(full=True)\n",
"V.set_outflux(-10.)\n",
"for i in range(10):\n",
" V.level = V.update_level(10.)\n",
" print(V.get_current_level())"
]
} }
], ],
"metadata": { "metadata": {

346
Untertweng_mit_Pegelregler.ipynb
View File

@@ -24,61 +24,58 @@
"source": [ "source": [
"# define constants\n", "# define constants\n",
"\n", "\n",
"#Turbine\n", " # for physics\n",
"Q_nenn = 0.85 # m³/s\n", "g = 9.81 # [m/s²] gravitational acceleration \n",
"p_nenn = pressure_conversion(10.6,'bar','Pa')\n", "rho = 1000. # [kg/m³] density of water \n",
"closing_time = 30. #s\n", "pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n",
"\n", "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n",
"# physics\n",
"g = 9.81 # gravitational acceleration [m/s²]\n",
"rho = 1000. # density of water [kg/m³]\n",
"\n", "\n",
"\n", "\n",
"# define controller constants\n", " # for Turbine\n",
"target_level = 8. # m\n", "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n",
"Kp = 0.1\n", "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n",
"Ti = 7.\n", "Tur_closingTime = 90. # [s] closing time of turbine\n",
"deadband_range = 0.05 # m\n",
"\n", "\n",
"\n", "\n",
"# pipeline\n", " # for PI controller\n",
"L = 535.+478. # length of pipeline [m]\n", "Con_targetLevel = 8. # [m]\n",
"D = 0.9 # pipe diameter [m]\n", "Con_K_p = 0.1 # [-] proportional constant of PI controller\n",
"A_pipe = D**2/4*np.pi # pipeline area\n", "Con_T_i = 1000. # [s] timespan in which a steady state error is corrected by the intergal term\n",
"h_pipe = 105 # hydraulic head without reservoir [m] \n", "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n",
"alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n",
"n = 50 # number of pipe segments in discretization\n",
"f_D = 0.014 # Darcy friction factor\n",
"c = 400. # 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 = 12000 # number of time steps after initial conditions\n",
"\n", "\n",
"\n",
" # for pipeline\n",
"Pip_length = (535.+478.) # [m] length of pipeline\n",
"Pip_dia = 0.9 # [m] diameter of pipeline\n",
"Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n",
"Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n",
"Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n",
"Pip_n_seg = 50 # [-] number of pipe segments in discretization\n",
"Pip_f_D = 0.014 # [-] Darcy friction factor\n",
"Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n",
" # derivatives of the pipeline constants\n", " # derivatives of the pipeline constants\n",
"dx = L/n # length of each pipe segment\n", "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n",
"dt = dx/c # timestep according to method of characterisitics\n", "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n",
"nn = n+1 # number of nodes\n", "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n",
"initial_level = target_level # water level in upstream reservoir [m]\n", "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n",
"pl_vec = np.arange(0,nn,1)*dx # pl = pipe-length. position of the nodes on the pipeline\n", "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n",
"t_vec = np.arange(0,nt+1)*dt # time vector\n",
"h_vec = np.arange(0,nn,1)*h_pipe/n # hydraulic head of pipeline at each node \n",
"\n", "\n",
"\n", "\n",
" # for reservoir\n",
"Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n",
"Res_area_out = Pip_area # [m²] outflux area of the reservoir, given by pipeline area\n",
"Res_level_crit_lo = 0. # [m] for yet-to-be-implemented warnings\n",
"Res_level_crit_hi = np.inf # [m] for yet-to-be-implemented warnings\n",
"Res_dt_approx = 1e-3 # [s] approx. timestep of reservoir time evolution to ensure numerical stability (see Res_nt why approx.)\n",
"Res_nt = max(1,int(Pip_dt//Res_dt_approx)) # [1] number of timesteps of the reservoir time evolution within one timestep of the pipeline\n",
"Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n",
"\n", "\n",
"# reservoir\n", " # for general simulation\n",
"# replace influx by vector\n", "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n",
"initial_flux = Q_nenn/1.1 # initial influx of volume to the reservoir [m³/s]\n", "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n",
"initial_pressure_unit = 'Pa' # DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", "simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n",
"conversion_pressure_unit = 'bar' # for pressure conversion in print statements and plot labels\n", "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n",
"area_base = 74. # total base are of the cuboid reservoir [m²] \n", "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n"
"area_outflux = A_pipe # outlfux area of the reservoir, given by pipeline area [m²]\n",
"critical_level_low = 0. # for yet-to-be-implemented warnings[m]\n",
"critical_level_high = np.inf # for yet-to-be-implemented warnings[m]\n",
"\n",
"# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n",
"nt_eRK4 = 100 # number of simulation steps of reservoir in between timesteps of pipeline \n",
"simulation_timestep = dt/nt_eRK4\n",
"\n",
"\n"
] ]
}, },
{ {
@@ -89,87 +86,101 @@
"source": [ "source": [
"# create objects\n", "# create objects\n",
"\n", "\n",
"V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,simulation_timestep)\n", "# Upstream reservoir\n",
"V.set_steady_state(initial_flux,initial_level,conversion_pressure_unit)\n", "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n",
"reservoir.set_steady_state(flux_init,level_init)\n",
"\n", "\n",
"# pipeline\n",
"pipe = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_n_seg,Pip_angle,Pip_f_D,Pip_pw_vel,Pip_dt,pUnit_conv,rho)\n",
"pipe.set_steady_state(flux_init,level_init,Res_area_base,Pip_x_vec,Pip_h_vec)\n",
"\n", "\n",
"pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n", "# downstream turbine\n",
"pipe.set_pressure_propagation_velocity(c)\n", "turbine = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime,Pip_dt,pUnit_conv)\n",
"pipe.set_number_of_timesteps(nt)\n", "turbine.set_steady_state(flux_init,pipe.get_current_pressure_distribution()[-1])\n",
"pipe.set_steady_state(initial_flux,initial_level,area_base,pl_vec,h_vec)\n",
"\n", "\n",
"initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n", "# influx setting turbine\n",
"turbine_in = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime,Pip_dt,pUnit_conv)\n",
"turbine_in.set_steady_state(flux_init,Tur_p_nenn)\n",
"\n", "\n",
"T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time,timestep=dt)\n", "# level controll\n",
"T1.set_steady_state(initial_flux,initial_pressure_turbine)\n", "level_control = PI_controller_class(Con_targetLevel,Con_deadbandRange,Con_K_p,Con_T_i,Pip_dt)\n",
"\n", "level_control.set_control_variable(turbine.get_current_LA(),display_warning=False)\n"
"T_in = Francis_Turbine(Q_nenn,p_nenn,closing_time/2,timestep=dt)\n",
"T_in.set_steady_state(initial_flux,p_nenn)\n",
"\n",
"Pegelregler = PI_controller_class(target_level,deadband_range,Kp,Ti,dt)\n",
"Pegelregler.control_variable = T1.get_current_LA()\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 11, "execution_count": 5,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# initialization for timeloop\n", "# initialization for Timeloop\n",
"\n", "\n",
"# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n",
"v_old = pipe.get_current_velocity_distribution()\n", "v_old = pipe.get_current_velocity_distribution()\n",
"v_min = pipe.get_current_velocity_distribution()\n",
"v_max = pipe.get_current_velocity_distribution()\n",
"Q_old = pipe.get_current_flux_distribution()\n",
"Q_min = pipe.get_current_flux_distribution()\n",
"Q_max = pipe.get_current_flux_distribution()\n",
"p_old = pipe.get_current_pressure_distribution()\n", "p_old = pipe.get_current_pressure_distribution()\n",
"p_min = pipe.get_current_pressure_distribution()\n",
"p_max = pipe.get_current_pressure_distribution()\n",
"\n",
"Q_in_vec = np.zeros_like(t_vec)\n",
"Q_in_vec[0] = flux_init\n",
"\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 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_res = np.zeros_like(t_vec)\n",
"v_boundary_tur = np.zeros_like(t_vec)\n", "v_boundary_tur = np.zeros_like(t_vec)\n",
"Q_boundary_res = np.zeros_like(t_vec)\n",
"Q_boundary_tur = np.zeros_like(t_vec)\n",
"p_boundary_res = np.zeros_like(t_vec)\n", "p_boundary_res = np.zeros_like(t_vec)\n",
"p_boundary_tur = np.zeros_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", "level_vec = np.full_like(t_vec,level_init) # level at the end of each pipeline timestep\n",
"level_vec = np.full(nt+1,initial_level) # level at the end of each pipeline timestep\n", "volume_vec = np.full_like(t_vec,reservoir.get_current_volume()) # volume at the end of each pipeline timestep\n",
"\n", "\n",
"# set the boundary conditions for the first timestep\n",
"v_boundary_res[0] = v_old[0]\n", "v_boundary_res[0] = v_old[0]\n",
"v_boundary_tur[0] = v_old[-1] \n", "v_boundary_tur[0] = v_old[-1] \n",
"Q_boundary_res[0] = Q_old[0]\n",
"Q_boundary_tur[0] = Q_old[-1]\n",
"p_boundary_res[0] = p_old[0]\n", "p_boundary_res[0] = p_old[0]\n",
"p_boundary_tur[0] = p_old[-1]\n", "p_boundary_tur[0] = p_old[-1]\n",
"\n", "\n",
"LA_soll_vec = np.full_like(t_vec,T1.get_current_LA())\n", "LA_soll_vec = np.full_like(t_vec,turbine.get_current_LA())\n",
"LA_ist_vec = np.full_like(t_vec,T1.get_current_LA())\n", "LA_ist_vec = np.full_like(t_vec,turbine.get_current_LA())\n",
"\n", "\n",
"LA_soll_vec2 = np.full_like(t_vec,T_in.get_current_LA())\n", "LA_soll_vec2 = np.full_like(t_vec,turbine_in.get_current_LA())\n",
"LA_soll_vec2[500:1000] = 0.\n", "LA_soll_vec2[500:] = 0\n",
"LA_soll_vec2[1000:1500] = 1. \n", "# LA_soll_vec2[500:1000] = 0.\n",
"LA_soll_vec2[1500:2000] = 0.\n", "# LA_soll_vec2[1000:1500] = 1. \n",
"LA_soll_vec2[2000:2500] = 0.5 \n" "# LA_soll_vec2[1500:2000] = 0.\n",
"# LA_soll_vec2[2000:2500] = 0.5 \n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 12, "execution_count": 6,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"%matplotlib qt5\n", "%matplotlib qt5\n",
"# time loop\n", "# Con_T_ime loop\n",
"\n", "\n",
"# create a figure and subplots to display the velocity and pressure distribution across the pipeline in each pipeline step\n", "# create a figure and subplots to display the velocity and pressure distribution across the pipeline in each pipeline step\n",
"fig1,axs1 = plt.subplots(2,1)\n", "fig1,axs1 = plt.subplots(2,1)\n",
"fig1.suptitle(str(0) +' s / '+str(round(t_vec[-1],2)) + ' s' )\n", "fig1.suptitle(str(0) +' s / '+str(round(t_vec[-1],2)) + ' s' )\n",
"axs1[0].set_title('Pressure distribution in pipeline')\n", "axs1[0].set_title('Pressure distribution in pipeline')\n",
"axs1[1].set_title('Velocity distribution in pipeline')\n",
"axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs1[0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", "axs1[0].set_ylabel(r'$p$ ['+pUnit_conv+']')\n",
"axs1[1].set_title('Flux distribution in pipeline')\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'$Q$ [$\\mathrm{m}^3 / \\mathrm{s}$]')\n",
"lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.')\n", "lo_p, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,pUnit_calc, pUnit_conv),marker='.')\n",
"lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n", "lo_q, = axs1[1].plot(Pip_x_vec,Q_old,marker='.')\n",
"lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n",
"lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n",
"lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n",
"lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n",
"\n",
"axs1[0].autoscale()\n", "axs1[0].autoscale()\n",
"axs1[1].autoscale()\n", "axs1[1].autoscale()\n",
"\n", "\n",
@@ -180,111 +191,150 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 13, "execution_count": 7,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# loop through time steps of the pipeline\n", "convergence_parameters = [p_old[-2],v_old[-2],Pip_dia,Pip_area,Pip_angle,Pip_f_D,Pip_pw_vel,rho,Pip_dt]\n",
"for it_pipe in range(1,pipe.nt+1):\n",
"\n", "\n",
" T_in.update_LA(LA_soll_vec2[it_pipe])\n", "# loop through Con_T_ime steps of the pipeline\n",
" T_in.set_pressure(p_nenn)\n", "for it_pipe in range(1,nt+1):\n",
" V.set_influx(T_in.get_current_Q())\n", "\n",
" turbine_in.update_LA(LA_soll_vec2[it_pipe])\n",
" turbine_in.set_pressure(Tur_p_nenn)\n",
" Q_in_vec[it_pipe] = turbine_in.get_current_Q()\n",
" reservoir.set_influx(Q_in_vec[it_pipe])\n",
"\n", "\n",
"# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\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", " # set initial condition for the reservoir Con_T_ime evolution calculted with e-RK4\n",
" V.set_pressure(p_old[0])\n", " reservoir.set_pressure(p_old[0],display_warning=False)\n",
" V.set_outflux(v_old[0]*area_outflux)\n", " reservoir.set_outflux(Q_old[0],display_warning=False)\n",
" # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n", " # calculate the Con_T_ime evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n",
" for it_res in range(nt_eRK4):\n", " for it_res in range(Res_nt):\n",
" V.timestep_reservoir_evolution() \n", " reservoir.timestep_reservoir_evolution() \n",
" level_vec[it_pipe] = V.get_current_level() \n", " level_vec[it_pipe] = reservoir.get_current_level() \n",
" volume_vec[it_pipe] = reservoir.get_current_volume() \n",
"\n", "\n",
" # get the control variable\n", " # get the control variable\n",
" Pegelregler.update_control_variable(level_vec[it_pipe])\n", " level_control.update_control_variable(level_vec[it_pipe])\n",
" LA_soll_vec[it_pipe] = Pegelregler.get_current_control_variable()\n", " LA_soll_vec[it_pipe] = level_control.get_current_control_variable()\n",
" \n", " \n",
" # change the Leitapparatöffnung based on the target value\n", " # change the Leitapparatöffnung based on the target value\n",
" T1.update_LA(LA_soll_vec[it_pipe])\n", " turbine.update_LA(LA_soll_vec[it_pipe])\n",
" LA_ist_vec[it_pipe] = T1.get_current_LA()\n", " LA_ist_vec[it_pipe] = turbine.get_current_LA()\n",
"\n", "\n",
" T1.set_pressure(p_old[-1])\n", " # set boundary condition for the next timestep of the characterisCon_T_ic method\n",
" # set boundary conditions for the next timestep of the characteristic method\n", " turbine.set_pressure(p_old[-1])\n",
" p_boundary_res[it_pipe] = V.get_current_pressure()\n", " convergence_parameters[0] = p_old[-2]\n",
" v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_current_Q()\n", " convergence_parameters[1] = v_old[-2]\n",
" turbine.converge(convergence_parameters)\n",
" p_boundary_res[it_pipe] = reservoir.get_current_pressure()\n",
" v_boundary_tur[it_pipe] = 1/Pip_area*turbine.get_current_Q()\n",
" Q_boundary_tur[it_pipe] = turbine.get_current_Q()\n",
"\n", "\n",
" # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n", " # the the boundary condition 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", " pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n",
" pipe.v[0] = (0.8*pipe.v[0]+0.2*reservoir.get_current_outflux()/Res_area_out)\n",
" p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\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", " v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n",
" Q_boundary_res[it_pipe] = pipe.get_current_flux_distribution()[0]\n",
"\n", "\n",
"\n", " # perform the next timestep via the characterisCon_T_ic method\n",
" # perform the next timestep via the characteristic method\n",
" pipe.timestep_characteristic_method()\n", " pipe.timestep_characteristic_method()\n",
"\n", "\n",
" # prepare for next loop\n", " # prepare for next loop\n",
" p_old = pipe.get_current_pressure_distribution()\n", " p_old = pipe.get_current_pressure_distribution()\n",
" v_old = pipe.get_current_velocity_distribution()\n", " v_old = pipe.get_current_velocity_distribution()\n",
" Q_old = pipe.get_current_flux_distribution()\n",
"\n",
"\n", "\n",
" # plot some stuff\n", " # plot some stuff\n",
" # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n", " # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n",
" lo_00.remove()\n", " lo_p.remove()\n",
" lo_01.remove()\n", " lo_pmin.remove()\n",
" # lo_02.remove()\n", " lo_pmax.remove()\n",
" lo_q.remove()\n",
" lo_qmin.remove()\n",
" lo_qmax.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(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.',c='blue')\n", " lo_p, = axs1[0].plot(Pip_x_vec,pipe.get_current_pressure_distribution(disp=True),marker='.',c='blue')\n",
" lo_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n", " lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n",
" # lo_02, = axs1[2].plot(level_vec_2,c='blue')\n", " lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n",
" lo_q, = axs1[1].plot(Pip_x_vec,pipe.get_current_flux_distribution(),marker='.',c='blue')\n",
" lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n",
" lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\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",
" fig1.canvas.draw()\n", " fig1.canvas.draw()\n",
" fig1.tight_layout()\n", " fig1.tight_layout()\n",
" fig1.show()\n", " fig1.show()\n",
" plt.pause(0.001) \n", " plt.pause(0.001) "
"\n",
" \n",
" "
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 14, "execution_count": 13,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# plot time evolution of boundary pressure and velocity as well as the reservoir level\n", "# plot Con_T_ime 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(1,1)\n",
"axs2[0,0].set_title('Pressure reservoir')\n", "axs2.set_title('Level and Volume reservoir')\n",
"axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,initial_pressure_unit, conversion_pressure_unit))\n", "axs2.plot(t_vec,level_vec,label='level')\n",
"axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[0,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", "axs2.set_ylabel(r'$h$ [m]')\n",
"x_twin_00 = axs2.twinx()\n",
"x_twin_00.set_ylabel(r'$V$ [$\\mathrm{m}^3$]')\n",
"x_twin_00.plot(t_vec,volume_vec)\n",
"axs2.legend()\n",
"\n", "\n",
"axs2[0,1].set_title('Velocity reservoir')\n", "fig2,axs2 = plt.subplots(1,1)\n",
"axs2[0,1].plot(t_vec,v_boundary_res)\n", "axs2.set_title('LA')\n",
"axs2[0,1].set_ylim(-2*Q_nenn,+2*Q_nenn)\n", "axs2.plot(t_vec,100*LA_soll_vec,label='Target')\n",
"axs2[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2.plot(t_vec,100*LA_ist_vec,label='Actual')\n",
"axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2.set_ylabel(r'$LA$ [%]')\n",
"axs2.legend()\n",
"\n", "\n",
"axs2[1,0].set_title('Pressure turbine')\n", "fig2,axs2 = plt.subplots(1,1)\n",
"axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,initial_pressure_unit, conversion_pressure_unit))\n", "axs2.set_title('Pressure reservoir and turbine')\n",
"axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2.plot(t_vec,pressure_conversion(p_boundary_res,pUnit_calc, pUnit_conv),label='Reservoir')\n",
"axs2[1,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", "axs2.plot(t_vec,pressure_conversion(p_boundary_tur,pUnit_calc, pUnit_conv),label='Turbine')\n",
"axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n",
"axs2.legend()\n",
"\n", "\n",
"axs2[1,1].set_title('Velocity turbine')\n", "fig2,axs2 = plt.subplots(1,1)\n",
"axs2[1,1].plot(t_vec,v_boundary_tur)\n", "axs2.set_title('Fluxes')\n",
"axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2.plot(t_vec,Q_boundary_res,label='Outflux')\n",
"axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", "axs2.plot(t_vec,Q_in_vec,label='Influx')\n",
"axs2.plot(t_vec,Q_boundary_tur,label='Flux Turbine')\n",
"axs2.set_ylim(-2*Tur_Q_nenn,+2*Tur_Q_nenn)\n",
"axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2.set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n",
"axs2.legend()\n",
"\n",
"fig2,axs2 = plt.subplots(1,1)\n",
"axs2.set_title('Min and Max Pressure')\n",
"axs2.plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n",
"axs2.plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n",
"axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n",
"\n",
"fig2,axs2 = plt.subplots(1,1)\n",
"axs2.set_title('Min and Max Fluxes')\n",
"axs2.plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n",
"axs2.plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n",
"axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs2.set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n",
"\n",
"# axs2[0,1].legend()\n",
"# axs2[1,0].legend()\n",
"# axs2[1,1].legend()\n",
"# # axs2[2,0].legend()\n",
"# # axs2[2,1].legend()\n",
"\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",
"\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", "fig2.tight_layout()\n",
"plt.show()" "plt.show()"
] ]