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

small changes for consistency, comments and a small fix in the convergence method of the turbine

Browse Source
This commit is contained in:
Brantegger Georg
2022-08-08 14:49:22 +02:00
parent 5a790d5ca5
commit 38c809ef49
10 changed files with 496 additions and 304 deletions
Download Patch File Download Diff File Expand all files Collapse all files

98
Ausgleichsbecken/Ausgleichsbecken_class_file.py
View File

@@ -14,7 +14,7 @@ def FODE_function(x_out,h,A,A_a,p,rho,g):
# based on the outflux formula by Andreas Malcherek # based on the outflux formula by Andreas Malcherek
# https://www.youtube.com/watch?v=8HO2LwqOhqQ # https://www.youtube.com/watch?v=8HO2LwqOhqQ
# adapted for a pressurized pipeline into which the reservoir effuses # adapted for a pressurized pipeline into which the reservoir effuses
# 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 ... Area_outflux # A_a ... Area_outflux
@@ -27,14 +27,14 @@ 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 display units are the same
# units are used to label graphs and disp 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$'
flux_unit = r'$\mathrm{m}^3/\mathrm{s}$' flux_unit = r'$\mathrm{m}^3/\mathrm{s}$'
level_unit = 'm' level_unit = 'm'
pressure_unit = 'Pa' pressure_unit = 'Pa' # DONT CHANGE needed for pressure conversion
time_unit = 's' time_unit = 's'
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$'
@@ -44,6 +44,7 @@ class Ausgleichsbecken_class:
density_unit_disp = 'kg/m³' density_unit_disp = 'kg/m³'
flux_unit_disp = 'm³/s' flux_unit_disp = 'm³/s'
level_unit_disp = 'm' level_unit_disp = 'm'
# pressure_unit_disp will be set within the __init__() method
time_unit_disp = 's' time_unit_disp = 's'
velocity_unit_disp = 'm/s' velocity_unit_disp = 'm/s'
volume_unit_disp = 'm³' volume_unit_disp = 'm³'
@@ -53,26 +54,37 @@ class Ausgleichsbecken_class:
# init # init
def __init__(self,area,area_outflux,timestep,pressure_unit_disp,level_min=0,level_max=np.inf,rho = 1000.): def __init__(self,area,area_outflux,timestep,pressure_unit_disp,level_min=0,level_max=np.inf,rho = 1000.):
"""
Creates a reservoir with given attributes in this order: \n
Base Area [m²] \n
Outflux Area [m²] \n
Simulation timestep [s] \n
Pressure unit for displaying [string] \n
Minimal level [m] \n
Maximal level [m] \n
Density of the liquid [kg/m³] \n
"""
#set initial attributes
self.area = area # base area of the cuboid reservoir self.area = area # base area of the cuboid reservoir
self.area_out = area_outflux # 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 - warning yet to be implemented
self.level_max = level_max # highest allowed water level self.level_max = level_max # highest allowed water level - warning yet to be implemented
self.pressure_unit_disp = pressure_unit_disp # pressure unit for displaying self.pressure_unit_disp = pressure_unit_disp # pressure unit for displaying
self.timestep = timestep # timestep in the time evolution method self.timestep = timestep # timestep in the time evolution method
# initialize for get_info # initialize for get_info() (if get_info() gets called before set_steady_state() is executed)
self.influx = "--" self.influx = -np.inf
self.outflux = "--" self.outflux = -np.inf
self.level = "--" self.level = -np.inf
self.pressure = "--" self.pressure = -np.inf
self.volume = "--" self.volume = -np.inf
# setter # setter - set attributes
def set_initial_level(self,initial_level): def set_initial_level(self,initial_level):
# sets the initial 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 == -np.inf:
self.level = initial_level self.level = initial_level
self.update_volume(set_flag=True) self.update_volume(set_flag=True)
else: else:
@@ -80,7 +92,7 @@ class Ausgleichsbecken_class:
def set_initial_pressure(self,initial_pressure): def set_initial_pressure(self,initial_pressure):
# sets the initial static pressure present at the outlet of the reservoir and should only be called during initialization # sets the initial static pressure present at the outlet of the reservoir and should only be called during initialization
if self.pressure == '--': if self.pressure == -np.inf:
self.pressure = initial_pressure self.pressure = initial_pressure
else: else:
raise Exception('Initial pressure was already set once. Use the .update_pressure(self) method to update pressure based current level.') raise Exception('Initial pressure was already set once. Use the .update_pressure(self) method to update pressure based current level.')
@@ -96,7 +108,7 @@ class Ausgleichsbecken_class:
if display_warning == True: if display_warning == True:
print('You are setting the outflux from the reservoir manually. \n \ print('You are setting the outflux from the reservoir manually. \n \
This is not an intended use of this method. \n \ This is not an intended use of this method. \n \
Refer to the timestep_reservoir_evolution() method instead.') Refer to the timestep_reservoir_evolution() or set_steady_state() method instead.')
self.outflux = outflux self.outflux = outflux
def set_level(self,level,display_warning=True): def set_level(self,level,display_warning=True):
@@ -104,7 +116,7 @@ class Ausgleichsbecken_class:
if display_warning == True: if display_warning == True:
print('You are setting the level of the reservoir manually. \n \ print('You are setting the level of the reservoir manually. \n \
This is not an intended use of this method. \n \ This is not an intended use of this method. \n \
Refer to the update_level() method instead.') Refer to the update_level() or set_steady_state() method instead.')
self.level = level self.level = level
def set_pressure(self,pressure,display_warning=True): def set_pressure(self,pressure,display_warning=True):
@@ -112,48 +124,53 @@ class Ausgleichsbecken_class:
if display_warning == True: if display_warning == True:
print('You are setting the pressure below the reservoir manually. \n \ print('You are setting the pressure below the reservoir manually. \n \
This is not an intended use of this method. \n \ This is not an intended use of this method. \n \
Refer to the update_pressure() method instead.') Refer to the update_pressure() or set_steady_state() method instead.')
self.pressure = pressure self.pressure = pressure
def set_volume(self,volume,display_warning=True): def set_volume(self,volume,display_warning=True):
# sets volume in reservoir
if display_warning == True: if display_warning == True:
print('You are setting the volume in the reservoir manually. \n \ print('You are setting the volume in the reservoir manually. \n \
This is not an intended use of this method. \n \ This is not an intended use of this method. \n \
Refer to the .update_volume() or set_initial_level() method instead.') Refer to the .update_volume() or set_initial_level() or set_steady_state() method instead.')
self.volume = volume self.volume = volume
def set_steady_state(self,ss_influx,ss_level): def set_steady_state(self,ss_influx,ss_level):
# set the steady state (ss) condition in which the net flux is zero # set the reservoir to 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_out) ss_outflux_vel = abs(ss_outflux/self.area_out)
# see confluence doc for explaination on how to arrive at the ss pressure formula
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)
# use setter methods to set the attributes to their steady state values
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) self.set_initial_pressure(ss_pressure)
self.set_outflux(ss_outflux,display_warning=False) self.set_outflux(ss_outflux,display_warning=False)
# getter # getter - return attributes
def get_info(self, full = False): def get_info(self, full = False):
# prints out the info on the current state of the reservoir
new_line = '\n' new_line = '\n'
p = pressure_conversion(self.pressure,self.pressure_unit,self.pressure_unit_disp) if self.pressure != np.inf:
outflux_vel = self.outflux/self.area_out p = pressure_conversion(self.pressure,self.pressure_unit,self.pressure_unit_disp)
if self.outflux != np.inf:
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_disp} {new_line}" f"Base area = {self.area:<10} {self.area_unit_disp} {new_line}"
f"Outflux area = {round(self.area_out,3):<10} {self.area_out_unit_disp} {new_line}" f"Outflux area = {round(self.area_out,3):<10} {self.area_outflux_unit_disp} {new_line}"
f"Current level = {self.level:<10} {self.level_unit_disp}{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_disp} {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_disp} {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_disp} {new_line}" f"Volume in reservoir = {self.volume:<10} {self.volume_unit_disp} {new_line}"
f"Current influx = {self.influx:<10} {self.flux_unit_disp} {new_line}" f"Current influx = {round(self.influx,3):<10} {self.flux_unit_disp} {new_line}"
f"Current outflux = {self.outflux:<10} {self.flux_unit_disp} {new_line}" f"Current outflux = {round(self.outflux,3):<10} {self.flux_unit_disp} {new_line}"
f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_disp} {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_disp} {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_disp} {new_line}" f"Simulation timestep = {self.timestep:<10} {self.time_unit_disp} {new_line}"
@@ -165,8 +182,8 @@ class Ausgleichsbecken_class:
f"----------------------------- {new_line}" f"----------------------------- {new_line}"
f"Current level = {self.level:<10} {self.level_unit_disp}{new_line}" f"Current level = {self.level:<10} {self.level_unit_disp}{new_line}"
f"Current volume = {self.volume:<10} {self.volume_unit_disp} {new_line}" f"Current volume = {self.volume:<10} {self.volume_unit_disp} {new_line}"
f"Current influx = {self.influx:<10} {self.flux_unit_disp} {new_line}" f"Current influx = {round(self.influx,3):<10} {self.flux_unit_disp} {new_line}"
f"Current outflux = {self.outflux:<10} {self.flux_unit_disp} {new_line}" f"Current outflux = {round(self.outflux,3):<10} {self.flux_unit_disp} {new_line}"
f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_disp} {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_disp} {new_line}" f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_disp} {new_line}"
f"----------------------------- {new_line}") f"----------------------------- {new_line}")
@@ -188,22 +205,26 @@ class Ausgleichsbecken_class:
def get_current_volume(self): def get_current_volume(self):
return self.volume return self.volume
# update methods # update methods - update attributes based on some parameter
def update_level(self,timestep,set_flag=False): def update_level(self,timestep,set_flag=False):
# 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
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
level_new = (self.level+delta_level) level_new = (self.level+delta_level)
# set flag is necessary because update_level() is used to get a halfstep value in the time evoultion
if set_flag == True: if set_flag == True:
self.set_level(level_new,display_warning=False) self.set_level(level_new,display_warning=False)
elif set_flag == False: elif set_flag == False:
return level_new return level_new
def update_pressure(self,set_flag=False): def update_pressure(self,set_flag=False):
influx_vel = abs(self.influx/self.area) # update pressure based on level and flux velocities
# see confluence doc for explaination
influx_vel = abs(self.influx/self.area)
outflux_vel = abs(self.outflux/self.area_out) 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)
# set flag for consistency with update_level()
if set_flag ==True: if set_flag ==True:
self.set_pressure(p_new,display_warning=False) self.set_pressure(p_new,display_warning=False)
elif set_flag == False: elif set_flag == False:
@@ -211,6 +232,7 @@ class Ausgleichsbecken_class:
def update_volume(self,set_flag=False): def update_volume(self,set_flag=False):
volume_new = self.level*self.area volume_new = self.level*self.area
# set flag for consistency with update_level()
if set_flag == True: if set_flag == True:
self.set_volume(volume_new,display_warning=False) self.set_volume(volume_new,display_warning=False)
elif set_flag == False: elif set_flag == False:
@@ -218,7 +240,9 @@ class Ausgleichsbecken_class:
#methods #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, level, pressure and volume based on current pipeline pressure and waterlevel in reservoir
# get some variables
dt = self.timestep dt = self.timestep
rho = self.density rho = self.density
g = self.g g = self.g
@@ -229,7 +253,8 @@ class Ausgleichsbecken_class:
h_hs = self.update_level(dt/2) h_hs = self.update_level(dt/2)
p = self.pressure p = self.pressure
p_hs = self.pressure + rho*g*(h_hs-h) p_hs = self.pressure + rho*g*(h_hs-h)
# explicit 4 step Runge Kutta
# perform explicit 4 step Runge Kutta
Y1 = yn Y1 = yn
Y2 = yn + dt/2*FODE_function(Y1,h,A,A_a,p,rho,g) Y2 = yn + dt/2*FODE_function(Y1,h,A,A_a,p,rho,g)
Y3 = yn + dt/2*FODE_function(Y2,h_hs,A,A_a,p_hs,rho,g) Y3 = yn + dt/2*FODE_function(Y2,h_hs,A,A_a,p_hs,rho,g)
@@ -237,6 +262,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))
# set/update the attributes to their new values
self.set_outflux(ynp1*A_a,display_warning=False) self.set_outflux(ynp1*A_a,display_warning=False)
self.update_level(dt,set_flag=True) self.update_level(dt,set_flag=True)
self.update_volume(set_flag=True) self.update_volume(set_flag=True)

155
Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb
View File

File diff suppressed because one or more lines are too long

108
Druckrohrleitung/Druckrohrleitung_class_file.py
View File

@@ -10,13 +10,15 @@ from functions.pressure_conversion import pressure_conversion
class Druckrohrleitung_class: class Druckrohrleitung_class:
# units # units
# make sure that units and display units are the same
# units are used to label graphs and disp units are used to have a bearable format when using pythons print()
acceleration_unit = r'$\mathrm{m}/\mathrm{s}^2$' acceleration_unit = r'$\mathrm{m}/\mathrm{s}^2$'
angle_unit = 'rad' angle_unit = 'rad'
area_unit = r'$\mathrm{m}^2$' area_unit = r'$\mathrm{m}^2$'
density_unit = r'$\mathrm{kg}/\mathrm{m}^3$' density_unit = r'$\mathrm{kg}/\mathrm{m}^3$'
flux_unit = r'$\mathrm{m}^3/\mathrm{s}$' flux_unit = r'$\mathrm{m}^3/\mathrm{s}$'
length_unit = 'm' length_unit = 'm'
pressure_unit = 'Pa' pressure_unit = 'Pa' # DONT CHANGE needed for pressure conversion
time_unit = 's' time_unit = 's'
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$'
@@ -27,33 +29,54 @@ class Druckrohrleitung_class:
density_unit_disp = 'kg/m³' density_unit_disp = 'kg/m³'
flux_unit_disp = 'm³/s' flux_unit_disp = 'm³/s'
length_unit_disp = 'm' length_unit_disp = 'm'
# pressure_unit_disp will be set within the __init__() method
time_unit_disp = 's' time_unit_disp = 's'
velocity_unit_disp = 'm/s' # for flux and pressure propagation velocity_unit_disp = 'm/s' # for flux and pressure propagation
volume_unit_disp = 'm³' volume_unit_disp = 'm³'
g = 9.81 g = 9.81 # m/s² gravitational acceleration
# init # init
def __init__(self,total_length,diameter,number_segments,pipeline_angle,Darcy_friction_factor,pw_vel,timestep,pressure_unit_disp,rho=1000): def __init__(self,total_length,diameter,pipeline_head,number_segments,Darcy_friction_factor,pw_vel,timestep,pressure_unit_disp,rho=1000):
self.length = total_length # total length of the pipeline """
self.dia = diameter # diameter of the pipeline Creates a reservoir with given attributes in this order: \n
self.n_seg = number_segments # number of segments for the method of characteristics Pipeline length [m] \n
self.angle = pipeline_angle # angle of the pipeline Pipeline diameter [m] \n
self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient Pipeline head [m] \n
self.c = pw_vel Number of pipeline segments [1] \n
Darcy friction factor [1] \n
Pressure wave velocity [m/s] \n
Simulation timestep [s] \n
Pressure unit for displaying [string] \n
Density of the liquid [kg/m³] \n
"""
self.length = total_length # total length of the pipeline
self.dia = diameter # diameter of the pipeline
self.head = pipeline_head # hydraulic head of the pipeline
self.n_seg = number_segments # number of segments for the method of characteristics
self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient
self.c = pw_vel # propagation velocity of pressure wave
self.dt = timestep self.dt = timestep
self.density = rho # density of the liquid in the pipeline self.density = rho # density of the liquid in the pipeline
self.A = (diameter/2)**2*np.pi # derivatives of input attributes
self.angle = np.arcsin(self.head/self.length) # angle of the pipeline
self.A = (diameter/2)**2*np.pi # crossectional area of the pipeline
self.dx = total_length/number_segments # length of each segment
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
self.h_vec = np.arange(0,(number_segments+1),1)*self.head/self.n_seg # vector giving the height difference from each node to the start of the pipeline
self.dx = total_length/number_segments # length of each segment self.pressure_unit_disp = pressure_unit_disp # pressure unit for displaying
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
self.pressure_unit_disp = pressure_unit_disp # setter - set attributes
def set_initial_pressure(self,pressure,display_warning=True):
# setter
def set_initial_pressure(self,pressure):
# initialize the pressure distribution in the pipeline # initialize the pressure distribution in the pipeline
if display_warning == True:
print('You are setting the pressure distribution in the pipeline manually. \n \
This is not an intended use of this method. \n \
Refer to the set_steady_state() method instead.')
# make sure the vector has the right size
if np.size(pressure) == 1: if np.size(pressure) == 1:
p0 = np.full_like(self.x_vec,pressure) p0 = np.full_like(self.x_vec,pressure)
elif np.size(pressure) == np.size(self.x_vec): elif np.size(pressure) == np.size(self.x_vec):
@@ -64,11 +87,18 @@ class Druckrohrleitung_class:
#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 = p0.copy() self.p_old = p0.copy()
self.p = p0.copy() self.p = p0.copy()
# initialize the vectors in which the minimal and maximal value of the pressure at each node are stores
self.p_min = p0.copy() self.p_min = p0.copy()
self.p_max = p0.copy() self.p_max = p0.copy()
def set_initial_flow_velocity(self,velocity): def set_initial_flow_velocity(self,velocity, display_warning=True):
# initialize the velocity distribution in the pipeline # initialize the velocity distribution in the pipeline
if display_warning == True:
print('You are setting the velocity distribution in the pipeline manually. \n \
This is not an intended use of this method. \n \
Refer to the set_steady_state() method instead.')
# make sure the vector has the right size
if np.size(velocity) == 1: if np.size(velocity) == 1:
v0 = np.full_like(self.x_vec,velocity) v0 = np.full_like(self.x_vec,velocity)
elif np.size(velocity) == np.size(self.x_vec): elif np.size(velocity) == np.size(self.x_vec):
@@ -79,6 +109,7 @@ class Druckrohrleitung_class:
#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 = v0.copy() self.v_old = v0.copy()
self.v = v0.copy() self.v = v0.copy()
# initialize the vectors in which the minimal and maximal value of the velocity at each node are stores
self.v_min = v0.copy() self.v_min = v0.copy()
self.v_max = v0.copy() self.v_max = v0.copy()
@@ -114,21 +145,19 @@ 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,x_vec,h_vec): def set_steady_state(self,ss_flux,ss_pressure_res):
# 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_like(self.x_vec,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_pressure = ss_pressure_res+(self.density*self.g*self.h_vec)-(self.f_D*self.x_vec/self.dia*self.density/2*ss_v0**2)
ss_v_out_res = abs(ss_flux/self.A)
ss_pressure_res = self.density*self.g*(ss_level_reservoir)+self.density*ss_v_out_res*(ss_v_in_res-ss_v_out_res)
ss_pressure = ss_pressure_res+(self.density*self.g*h_vec)-(self.f_D*x_vec/self.dia*self.density/2*ss_v0**2)
self.set_initial_flow_velocity(ss_v0) # set the initial conditions
self.set_initial_pressure(ss_pressure) self.set_initial_flow_velocity(ss_v0,display_warning=False)
self.set_initial_pressure(ss_pressure,display_warning=False)
# getter # getter - return attributes
def get_info(self): def get_info(self):
new_line = '\n' new_line = '\n'
angle_deg = round(self.angle/np.pi*180,3) angle_deg = round(self.angle/np.pi*180,3)
@@ -139,6 +168,7 @@ class Druckrohrleitung_class:
f"----------------------------- {new_line}" f"----------------------------- {new_line}"
f"Length = {self.length:<10} {self.length_unit_disp} {new_line}" f"Length = {self.length:<10} {self.length_unit_disp} {new_line}"
f"Diameter = {self.dia:<10} {self.length_unit_disp} {new_line}" f"Diameter = {self.dia:<10} {self.length_unit_disp} {new_line}"
f"Hydraulic head = {self.head:<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_disp} {new_line}" f"Length per segments = {self.dx:<10} {self.length_unit_disp} {new_line}"
@@ -148,17 +178,16 @@ class Druckrohrleitung_class:
f"Density of liquid = {self.density:<10} {self.density_unit_disp} {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_disp} {new_line}" f"Pressure wave vel. = {self.c:<10} {self.velocity_unit_disp} {new_line}"
f"Simulation timestep = {self.dt:<10} {self.time_unit_disp} {new_line}" f"Simulation timestep = {self.dt:<10} {self.time_unit_disp} {new_line}"
f"Number of timesteps = {self.nt:<10} {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,disp=False): def get_current_pressure_distribution(self,disp_flag=False):
if disp == True: # disp_flag if one wants to directly plot the return of this method
if disp_flag == True: # convert to pressure unit disp
return pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_disp) return pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_disp)
elif disp == False: elif disp_flag == False: # stay in Pa
return self.p return self.p
def get_current_velocity_distribution(self): def get_current_velocity_distribution(self):
@@ -167,16 +196,16 @@ class Druckrohrleitung_class:
def get_current_flux_distribution(self): def get_current_flux_distribution(self):
return self.v*self.A return self.v*self.A
def get_lowest_pressure_per_node(self,disp=False): def get_lowest_pressure_per_node(self,disp_flag=False):
if disp == True: if disp_flag == True: # convert to pressure unit disp
return pressure_conversion(self.p_min,self.pressure_unit,self.pressure_unit_disp) return pressure_conversion(self.p_min,self.pressure_unit,self.pressure_unit_disp)
elif disp == False: elif disp_flag == False: # stay in Pa
return self.p_min return self.p_min
def get_highest_pressure_per_node(self,disp=False): def get_highest_pressure_per_node(self,disp_flag=False):
if disp == True: if disp_flag == True: # convert to pressure unit disp
return pressure_conversion(self.p_max,self.pressure_unit,self.pressure_unit_disp) return pressure_conversion(self.p_max,self.pressure_unit,self.pressure_unit_disp)
elif disp == False: elif disp_flag == False: # stay in Pa
return self.p_max return self.p_max
def get_lowest_velocity_per_node(self): def get_lowest_velocity_per_node(self):
@@ -194,14 +223,15 @@ class Druckrohrleitung_class:
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
# constants for cleaner formula
nn = self.n_seg+1 # number of nodes nn = self.n_seg+1 # number of nodes
rho = self.density # density of liquid rho = self.density # density of liquid
c = self.c # pressure propagation velocity c = self.c # pressure propagation velocity
f_D = self.f_D # Darcy friction coefficient f_D = self.f_D # Darcy friction coefficient
dt = self.dt # timestep dt = self.dt # timestep
D = self.dia # pipeline diametet D = self.dia # pipeline diameter
g = self.g # graviational acceleration g = self.g # graviational acceleration
alpha = self.angle # pipeline angle alpha = self.angle # pipeline angle

151
Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb
View File

@@ -34,36 +34,32 @@
"pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \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", "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n",
"\n", "\n",
"\n",
" # for Turbine\n", " # for Turbine\n",
"Tur_Q_nenn = 0.85 # [m³/s] nominal flux of 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_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n",
"Tur_closingTime = 90. # [s] closing time of turbine\n", "Tur_closingTime = 90. # [s] closing time of turbine\n",
"\n",
"\n", "\n",
" # for PI controller\n", " # for PI controller\n",
"Con_targetLevel = 8. # [m]\n", "Con_targetLevel = 8. # [m]\n",
"Con_K_p = 0.1 # [-] proportional constant of PI controller\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_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", "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n",
"\n", "\n",
"\n",
" # for pipeline\n", " # for pipeline\n",
"Pip_length = (535.+478.) # [m] length of pipeline\n", "Pip_length = (535.+478.) # [m] length of pipeline\n",
"Pip_dia = 0.9 # [m] diameter 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_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_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_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_n_seg = 50 # [-] number of pipe segments in discretization\n",
"Pip_f_D = 0.014 # [-] Darcy friction factor\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", "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",
"Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\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_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_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_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", "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", "\n",
" # for reservoir\n", " # for reservoir\n",
"Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n", "Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n",
@@ -75,38 +71,68 @@
"Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n", "Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n",
"\n", "\n",
" # for general simulation\n", " # for general simulation\n",
"flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \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", "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", "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", "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" "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": 3, "execution_count": 3,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The pipeline has the following attributes: \n",
"----------------------------- \n",
"Length = 1013.0 m \n",
"Diameter = 0.9 m \n",
"Hydraulic head = 105.0 m \n",
"Number of segments = 50 \n",
"Number of nodes = 51 \n",
"Length per segments = 20.26 m \n",
"Pipeline angle = 0.104 rad \n",
"Pipeline angle = 5.95° \n",
"Darcy friction factor = 0.014 \n",
"Density of liquid = 1000.0 kg/m³ \n",
"Pressure wave vel. = 500.0 m/s \n",
"Simulation timestep = 0.04052 s \n",
"----------------------------- \n",
"Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object\n",
"The pipeline has the following attributes: \n",
"----------------------------- \n",
"Length = 1013.0 m \n",
"Diameter = 0.9 m \n",
"Hydraulic head = 105.0 m \n",
"Number of segments = 50 \n",
"Number of nodes = 51 \n",
"Length per segments = 20.26 m \n",
"Pipeline angle = 0.104 rad \n",
"Pipeline angle = 5.95° \n",
"Darcy friction factor = 0.014 \n",
"Density of liquid = 1000.0 kg/m³ \n",
"Pressure wave vel. = 500.0 m/s \n",
"Simulation timestep = 0.04052 s \n",
"----------------------------- \n",
"Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object\n"
]
}
],
"source": [ "source": [
"# create objects\n", "# create objects\n",
"\n", "\n",
"# Upstream reservoir\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 = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,pUnit_conv,Res_level_crit_lo,Res_level_crit_hi,rho)\n",
"reservoir.set_steady_state(flux_init,level_init)\n", "reservoir.set_steady_state(flux_init,level_init)\n",
"\n", "\n",
"# pipeline\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 = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_head,Pip_n_seg,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" "pipe.set_steady_state(flux_init,reservoir.get_current_pressure())\n",
] "pipe.get_info()\n"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"reservoir.get_info(full=True)\n",
"pipe.get_info(full=True)"
] ]
}, },
{ {
@@ -169,7 +195,46 @@
"cell_type": "code", "cell_type": "code",
"execution_count": 6, "execution_count": 6,
"metadata": {}, "metadata": {},
"outputs": [], "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.0 m\n",
"Critical level low = 0.0 m \n",
"Critical level high = inf m \n",
"Volume in reservoir = 592.0 m³ \n",
"Current influx = 0.773 m³/s \n",
"Current outflux = 0.773 m³/s \n",
"Current outflux vel = 1.215 m/s \n",
"Current pipe pressure = 7.854 mWS \n",
"Simulation timestep = 0.001013 s \n",
"Density of liquid = 1000.0 kg/m³ \n",
"----------------------------- \n",
"\n",
"The pipeline has the following attributes: \n",
"----------------------------- \n",
"Length = 1013.0 m \n",
"Diameter = 0.9 m \n",
"Hydraulic head = 105.0 m \n",
"Number of segments = 50 \n",
"Number of nodes = 51 \n",
"Length per segments = 20.26 m \n",
"Pipeline angle = 0.104 rad \n",
"Pipeline angle = 5.95° \n",
"Darcy friction factor = 0.014 \n",
"Density of liquid = 1000.0 kg/m³ \n",
"Pressure wave vel. = 500.0 m/s \n",
"Simulation timestep = 0.04052 s \n",
"----------------------------- \n",
"Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object\n"
]
}
],
"source": [ "source": [
"for it_pipe in range(1,nt+1):\n", "for it_pipe in range(1,nt+1):\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",
@@ -213,12 +278,12 @@
" plt.pause(0.000001)\n", " plt.pause(0.000001)\n",
"\n", "\n",
"reservoir.get_info(full=True)\n", "reservoir.get_info(full=True)\n",
"pipe.get_info(full=True)" "pipe.get_info()"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 12, "execution_count": 7,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [

14
Regler/Pegelregler_test.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 27, "execution_count": 1,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -23,7 +23,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 28, "execution_count": 2,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -85,7 +85,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 29, "execution_count": 3,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -93,7 +93,7 @@
"offset_pressure = pressure_conversion(Pip_head,'mws',pUnit_calc)\n", "offset_pressure = pressure_conversion(Pip_head,'mws',pUnit_calc)\n",
"\n", "\n",
"# Upstream reservoir\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 = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,pUnit_conv,Res_level_crit_lo,Res_level_crit_hi,rho)\n",
"reservoir.set_steady_state(flux_init,level_init)\n", "reservoir.set_steady_state(flux_init,level_init)\n",
"\n", "\n",
"# downstream turbine\n", "# downstream turbine\n",
@@ -108,7 +108,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 30, "execution_count": 4,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -124,7 +124,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 31, "execution_count": 5,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -181,7 +181,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 32, "execution_count": 6,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [

23
Regler/Regler_class_file.py
View File

@@ -1,6 +1,7 @@
import numpy as np import numpy as np
#based on https://en.wikipedia.org/wiki/PID_controller#Discrete_implementation #based on https://en.wikipedia.org/wiki/PID_controller#Discrete_implementation
# performance parameters for controllers
def trap_int(vec,timestep): def trap_int(vec,timestep):
# numerical integration via the trapeziod rule to calculate the performance parameters # numerical integration via the trapeziod rule to calculate the performance parameters
l = np.size(vec) l = np.size(vec)
@@ -41,6 +42,7 @@ def ITAE_fun(error_history,timestep):
itae = trap_int(np.abs(e),dt) itae = trap_int(np.abs(e),dt)
return itae return itae
# P controller
class P_controller_class: class P_controller_class:
# def __init__(self,setpoint,proportionality_constant): # def __init__(self,setpoint,proportionality_constant):
# self.SP = setpoint # self.SP = setpoint
@@ -72,14 +74,14 @@ class P_controller_class:
def __init__(self): def __init__(self):
pass pass
# PI controller
class PI_controller_class: class PI_controller_class:
# init # init
def __init__(self,setpoint,deadband,proportionality_constant,Ti,timestep,lower_limit=0.,upper_limit=1.): def __init__(self,setpoint,deadband,proportionality_constant,Ti,timestep,lower_limit=0.,upper_limit=1.):
self.SP = setpoint self.SP = setpoint
self.db = deadband self.db = deadband
self.Kp = proportionality_constant self.Kp = proportionality_constant
self.Ti = Ti # integration time self.Ti = Ti # ~integration time
self.dt = timestep self.dt = timestep
# 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]
@@ -88,15 +90,13 @@ class PI_controller_class:
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: 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') Consider using the .update_controll_variable() method instead.')
self.control_variable = control_variable self.control_variable = control_variable
# getter # getter
@@ -167,15 +167,14 @@ class PI_controller_class:
# only if that is the case, change control variable # only if that is the case, change control variable
if abs(self.error) > self.db: if abs(self.error) > self.db:
new_control = cv+Kp*(e0-e1)+dt/Ti*e0 new_control = cv+Kp*(e0-e1)+dt/Ti*e0
# ensure that the controll variable stays within the predefined limits
if new_control < self.cv_lower_limit:
new_control = self.cv_lower_limit
if new_control > self.cv_upper_limit:
new_control = self.cv_upper_limit
else: else:
new_control = cv new_control = cv
# ensure that the controll variable stays within the predefined limits
if new_control < self.cv_lower_limit:
new_control = self.cv_lower_limit
if new_control > self.cv_upper_limit:
new_control = self.cv_upper_limit
# set the control variable attribute # set the control variable attribute
self.set_control_variable(new_control,display_warning=False) self.set_control_variable(new_control,display_warning=False)

48
Turbinen/Turbinen_class_file.py
View File

@@ -11,8 +11,8 @@ sys.path.append(parent)
from functions.pressure_conversion import pressure_conversion 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 display units are the same
# units are used to label graphs and disp 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}$'
@@ -32,7 +32,16 @@ class Francis_Turbine:
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,timestep,pressure_unit_disp): def __init__(self,Q_nenn,p_nenn,t_closing,timestep,pressure_unit_disp):
"""
Creates a turbine with given attributes in this order: \n
Nominal flux [m³/s] \n
Nominal pressure [Pa] \n
Closing time [s] \n
Simulation timestep [s] \n
Pressure unit for displaying [string] \n
"""
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
@@ -42,21 +51,21 @@ class Francis_Turbine:
self.pressure_unit_disp = pressure_unit_disp 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()
self.p = pressure_conversion(-1,self.pressure_unit_disp,self.pressure_unit) self.p = -np.inf
self.Q = -1. self.Q = -np.inf
self.LA = -0.01 self.LA = -np.inf
# setter # setter - set attributes
def set_LA(self,LA,display_warning=True): 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 # warn user, that the .set_LA() method should not be used ot set LA manually
if display_warning == True: if display_warning == True:
print('You are setting the guide vane opening of the turbine manually. \n \ print('You are setting the guide vane opening of the turbine manually. \n \
This is not an intended use of this method. \n \ This is not an intended use of this method. \n \
Refer to the .update_LA() method instead.') Refer to the .update_LA() method instead.')
# set Leitapparatöffnung
self.LA = LA
def set_pressure(self,pressure): def set_pressure(self,pressure):
# set pressure in front of the turbine # set pressure in front of the turbine
@@ -70,7 +79,7 @@ class Francis_Turbine:
raise Exception('LA out of range [0;1]') raise Exception('LA out of range [0;1]')
self.set_LA(ss_LA,display_warning=False) self.set_LA(ss_LA,display_warning=False)
#getter #getter - get attributes
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
# of the turbine and the Leitapparatöffnung # of the turbine and the Leitapparatöffnung
@@ -83,8 +92,11 @@ class Francis_Turbine:
def get_current_LA(self): def get_current_LA(self):
return self.LA return self.LA
def get_current_pressure(self): def get_current_pressure(self,disp_flag=True):
return pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_disp) if disp_flag == True:
return pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_disp)
else:
return self.p
def get_info(self, full = False): def get_info(self, full = False):
new_line = '\n' new_line = '\n'
@@ -135,7 +147,9 @@ class Francis_Turbine:
# methods # methods
def converge(self,convergence_parameters): def converge(self,convergence_parameters):
# 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 # 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
# new pressure from the forward characteristic are not compatible. # new pressure from the forward characteristic are not perfectly compatible.
# Therefore, iterate the flux and the pressure so long, until they converge
eps = 1e-12 # convergence criterion: iteration change < eps eps = 1e-12 # convergence criterion: iteration change < eps
iteration_change = 1. # change in Q from one iteration to the next iteration_change = 1. # change in Q from one iteration to the next
i = 0 # safety variable. break loop if it exceeds 1e6 iterations i = 0 # safety variable. break loop if it exceeds 1e6 iterations
@@ -150,20 +164,18 @@ class Francis_Turbine:
rho = convergence_parameters[7] # density of the liquid rho = convergence_parameters[7] # density of the liquid
dt = convergence_parameters[8] # timestep of the characteristic method dt = convergence_parameters[8] # timestep of the characteristic method
p_old = self.get_current_pressure()
Q_old = self.get_current_Q() Q_old = self.get_current_Q()
v_old = Q_old/area_pipe v_old = Q_old/area_pipe
while iteration_change > eps: while iteration_change > eps:
self.set_pressure(p_old) 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
self.set_pressure(p_new)
Q_new = self.get_current_Q() Q_new = self.get_current_Q()
v_new = Q_new/area_pipe 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) iteration_change = abs(Q_old-Q_new)
Q_old = Q_new.copy() Q_old = Q_new.copy()
p_old = p_new.copy()
v_old = v_new.copy() v_old = v_new.copy()
i = i+1 i = i+1
if i == 1e6: if i == 1e6:

116
Turbinen/Turbinen_test_steady_state.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 8, "execution_count": null,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -23,7 +23,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 9, "execution_count": null,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -35,12 +35,10 @@
"pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \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", "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n",
"\n", "\n",
"\n",
" # for Turbine\n", " # for Turbine\n",
"Tur_Q_nenn = 0.85 # [m³/s] nominal flux of 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_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n",
"Tur_closingTime = 90. # [s] closing time of turbine\n", "Tur_closingTime = 90. # [s] closing time of turbine\n",
"\n",
"\n", "\n",
" # for PI controller\n", " # for PI controller\n",
"Con_targetLevel = 8. # [m]\n", "Con_targetLevel = 8. # [m]\n",
@@ -48,23 +46,21 @@
"Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\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", "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n",
"\n", "\n",
"\n",
" # for pipeline\n", " # for pipeline\n",
"Pip_length = (535.+478.) # [m] length of pipeline\n", "Pip_length = (535.+478.) # [m] length of pipeline\n",
"Pip_dia = 0.9 # [m] diameter 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_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_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_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_n_seg = 50 # [-] number of pipe segments in discretization\n",
"Pip_f_D = 0.014 # [-] Darcy friction factor\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", "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",
"Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\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_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_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_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", "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", "\n",
" # for reservoir\n", " # for reservoir\n",
"Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n", "Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n",
@@ -76,16 +72,16 @@
"Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n", "Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n",
"\n", "\n",
" # for general simulation\n", " # for general simulation\n",
"flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \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", "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", "simTime_target = 100. # [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", "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" "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": 10, "execution_count": null,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -96,8 +92,8 @@
"reservoir.set_steady_state(flux_init,level_init)\n", "reservoir.set_steady_state(flux_init,level_init)\n",
"\n", "\n",
"# pipeline\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 = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_head,Pip_n_seg,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", "pipe.set_steady_state(flux_init,reservoir.get_current_pressure())\n",
"\n", "\n",
"# downstream turbine\n", "# downstream turbine\n",
"turbine = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime,Pip_dt,pUnit_conv)\n", "turbine = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime,Pip_dt,pUnit_conv)\n",
@@ -114,7 +110,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 11, "execution_count": null,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -133,22 +129,22 @@
"Q_in_vec = np.zeros_like(t_vec)\n", "Q_in_vec = np.zeros_like(t_vec)\n",
"Q_in_vec[0] = flux_init\n", "Q_in_vec[0] = flux_init\n",
"\n", "\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_res = np.zeros_like(t_vec)\n",
"Q_boundary_tur = 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",
"level_vec = np.full_like(t_vec,level_init) # level at the end of each pipeline timestep\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", "volume_vec = np.full_like(t_vec,reservoir.get_current_volume()) # volume at the end of each pipeline timestep\n",
"\n", "\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_res[0] = Q_old[0]\n",
"Q_boundary_tur[0] = Q_old[-1]\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,turbine.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,turbine.get_current_LA())\n", "LA_ist_vec = np.full_like(t_vec,turbine.get_current_LA())\n",
@@ -158,7 +154,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 12, "execution_count": null,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -176,8 +172,8 @@
"axs1[1].set_ylabel(r'$Q$ [$\\mathrm{m}^3 / \\mathrm{s}$]')\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_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_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_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp_flag=True),c='red')\n",
"lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", "lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp_flag=True),c='red')\n",
"lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),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", "lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n",
"\n", "\n",
@@ -191,7 +187,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 13, "execution_count": null,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -212,7 +208,7 @@
" # calculate the Con_T_ime 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(Res_nt):\n", " for it_res in range(Res_nt):\n",
" reservoir.timestep_reservoir_evolution() \n", " reservoir.timestep_reservoir_evolution() \n",
" level_vec[it_pipe] = reservoir.get_current_level() \n", " level_vec[it_pipe] = reservoir.get_current_level() \n",
" volume_vec[it_pipe] = reservoir.get_current_volume() \n", " volume_vec[it_pipe] = reservoir.get_current_volume() \n",
"\n", "\n",
" # get the control variable\n", " # get the control variable\n",
@@ -247,7 +243,6 @@
" v_old = pipe.get_current_velocity_distribution()\n", " v_old = pipe.get_current_velocity_distribution()\n",
" Q_old = pipe.get_current_flux_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_p.remove()\n", " lo_p.remove()\n",
@@ -257,9 +252,9 @@
" lo_qmin.remove()\n", " lo_qmin.remove()\n",
" lo_qmax.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_p, = axs1[0].plot(Pip_x_vec,pipe.get_current_pressure_distribution(disp=True),marker='.',c='blue')\n", " lo_p, = axs1[0].plot(Pip_x_vec,pipe.get_current_pressure_distribution(disp_flag=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_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp_flag=True),c='red')\n",
" lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", " lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp_flag=True),c='red')\n",
" lo_q, = axs1[1].plot(Pip_x_vec,pipe.get_current_flux_distribution(),marker='.',c='blue')\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_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", " lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n",
@@ -272,7 +267,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 14, "execution_count": null,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -316,8 +311,8 @@
"\n", "\n",
"fig2,axs2 = plt.subplots(1,1)\n", "fig2,axs2 = plt.subplots(1,1)\n",
"axs2.set_title('Min and Max Pressure')\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_lowest_pressure_per_node(disp_flag=True),c='red')\n",
"axs2.plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", "axs2.plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp_flag=True),c='red')\n",
"axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n", "axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n",
"\n", "\n",
@@ -328,13 +323,6 @@
"axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs2.set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n", "axs2.set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n",
"\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", "fig2.tight_layout()\n",
"plt.show()" "plt.show()"
] ]

81
Untertweng_mit_Pegelregler.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 1, "execution_count": 4,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -18,7 +18,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 2, "execution_count": 5,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -30,12 +30,10 @@
"pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \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", "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n",
"\n", "\n",
"\n",
" # for Turbine\n", " # for Turbine\n",
"Tur_Q_nenn = 0.85 # [m³/s] nominal flux of 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_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n",
"Tur_closingTime = 90. # [s] closing time of turbine\n", "Tur_closingTime = 90. # [s] closing time of turbine\n",
"\n",
"\n", "\n",
" # for PI controller\n", " # for PI controller\n",
"Con_targetLevel = 8. # [m]\n", "Con_targetLevel = 8. # [m]\n",
@@ -43,23 +41,21 @@
"Con_T_i = 1000. # [s] timespan in which a steady state error is corrected by the intergal term\n", "Con_T_i = 1000. # [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", "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n",
"\n", "\n",
"\n",
" # for pipeline\n", " # for pipeline\n",
"Pip_length = (535.+478.) # [m] length of pipeline\n", "Pip_length = (535.+478.) # [m] length of pipeline\n",
"Pip_dia = 0.9 # [m] diameter 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_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_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_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_n_seg = 50 # [-] number of pipe segments in discretization\n",
"Pip_f_D = 0.014 # [-] Darcy friction factor\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", "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",
"Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\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_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_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_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", "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", "\n",
" # for reservoir\n", " # for reservoir\n",
"Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n", "Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n",
@@ -71,28 +67,28 @@
"Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n", "Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n",
"\n", "\n",
" # for general simulation\n", " # for general simulation\n",
"flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \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", "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", "simTime_target = 100. # [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", "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" "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": 3, "execution_count": 6,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# create objects\n", "# create objects\n",
"\n", "\n",
"# Upstream reservoir\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 = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,pUnit_conv,Res_level_crit_lo,Res_level_crit_hi,rho)\n",
"reservoir.set_steady_state(flux_init,level_init)\n", "reservoir.set_steady_state(flux_init,level_init)\n",
"\n", "\n",
"# pipeline\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 = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_head,Pip_n_seg,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", "pipe.set_steady_state(flux_init,reservoir.get_current_pressure())\n",
"\n", "\n",
"# downstream turbine\n", "# downstream turbine\n",
"turbine = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime,Pip_dt,pUnit_conv)\n", "turbine = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime,Pip_dt,pUnit_conv)\n",
@@ -109,7 +105,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 5, "execution_count": 7,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -158,7 +154,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 6, "execution_count": 9,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -176,8 +172,8 @@
"axs1[1].set_ylabel(r'$Q$ [$\\mathrm{m}^3 / \\mathrm{s}$]')\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_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_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_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp_flag=True),c='red')\n",
"lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", "lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp_flag=True),c='red')\n",
"lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),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", "lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n",
"\n", "\n",
@@ -191,7 +187,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 7, "execution_count": 10,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -247,7 +243,6 @@
" v_old = pipe.get_current_velocity_distribution()\n", " v_old = pipe.get_current_velocity_distribution()\n",
" Q_old = pipe.get_current_flux_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_p.remove()\n", " lo_p.remove()\n",
@@ -257,9 +252,9 @@
" lo_qmin.remove()\n", " lo_qmin.remove()\n",
" lo_qmax.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_p, = axs1[0].plot(Pip_x_vec,pipe.get_current_pressure_distribution(disp=True),marker='.',c='blue')\n", " lo_p, = axs1[0].plot(Pip_x_vec,pipe.get_current_pressure_distribution(disp_flag=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_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp_flag=True),c='red')\n",
" lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", " lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp_flag=True),c='red')\n",
" lo_q, = axs1[1].plot(Pip_x_vec,pipe.get_current_flux_distribution(),marker='.',c='blue')\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_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", " lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n",
@@ -272,7 +267,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 13, "execution_count": 11,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -316,8 +311,8 @@
"\n", "\n",
"fig2,axs2 = plt.subplots(1,1)\n", "fig2,axs2 = plt.subplots(1,1)\n",
"axs2.set_title('Min and Max Pressure')\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_lowest_pressure_per_node(disp_flag=True),c='red')\n",
"axs2.plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", "axs2.plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp_flag=True),c='red')\n",
"axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n", "axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n",
"\n", "\n",

2
functions/pressure_conversion.py
View File

@@ -27,8 +27,6 @@ def pa_to_torr(p):
def pa_to_atm(p): def pa_to_atm(p):
return p*1/(101.325*1e3) return p*1/(101.325*1e3)
# converstion function
def pa_to_psi(p): def pa_to_psi(p):
return p/6894.8 return p/6894.8