diff --git a/Ausgleichsbecken/Ausgleichsbecken_class_file.py b/Ausgleichsbecken/Ausgleichsbecken_class_file.py index 970f7fa..31af119 100644 --- a/Ausgleichsbecken/Ausgleichsbecken_class_file.py +++ b/Ausgleichsbecken/Ausgleichsbecken_class_file.py @@ -1,6 +1,8 @@ -import os -import sys -from logging import exception +# import modules for general use +import os # to import functions from other folders +import sys # to import functions from other folders +from logging import \ + exception # to throw an exception when a specific condition is met import numpy as np @@ -13,10 +15,12 @@ from functions.pressure_conversion import pressure_conversion def FODE_function(x_out,h,A,A_a,p,rho,g): # (FODE ... first order differential equation) + # describes the change in outflux velocity from a reservoir # based on the outflux formula by Andreas Malcherek # https://www.youtube.com/watch?v=8HO2LwqOhqQ # adapted for a pressurized pipeline into which the reservoir effuses - # and flow direction + # and flow direction + # see documentation in word-file # x_out ... effusion velocity # h ... level in the reservoir # A_a ... Area_outflux @@ -29,14 +33,14 @@ def FODE_function(x_out,h,A,A_a,p,rho,g): class Ausgleichsbecken_class: # 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() + # make sure that units and display units are the same! + # units are used to label graphs and disp units are used to have good formatting when using pythons print() area_unit = r'$\mathrm{m}^2$' area_outflux_unit = r'$\mathrm{m}^2$' density_unit = r'$\mathrm{kg}/\mathrm{m}^3$' flux_unit = r'$\mathrm{m}^3/\mathrm{s}$' level_unit = 'm' - pressure_unit = 'Pa' # DONT CHANGE needed for pressure conversion + pressure_unit = 'Pa' # !DO NOT CHANGE! needed for pressure conversion time_unit = 's' velocity_unit = r'$\mathrm{m}/\mathrm{s}$' volume_unit = r'$\mathrm{m}^3$' @@ -55,6 +59,7 @@ class Ausgleichsbecken_class: # init + # see docstring below 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 @@ -62,8 +67,8 @@ class Ausgleichsbecken_class: Outflux Area [m²] \n Simulation timestep [s] \n Pressure unit for displaying [string] \n - Minimal level [m] \n - Maximal level [m] \n + Minimum level [m] \n + Maximum level [m] \n Density of the liquid [kg/m³] \n """ #set initial attributes @@ -75,7 +80,8 @@ class Ausgleichsbecken_class: self.pressure_unit_disp = pressure_unit_disp # pressure unit for displaying self.timestep = timestep # timestep in the time evolution method - # initialize for get_info() (if get_info() gets called before set_steady_state() is executed) + # initialize for get_info() (if get_info() gets called before set_steady_state() was ever executed) + # is also used to check if set_steady_state() was ever executed self.influx = -np.inf self.outflux = -np.inf self.level = -np.inf @@ -97,7 +103,7 @@ class Ausgleichsbecken_class: if self.pressure == -np.inf: 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.') + raise Exception('Initial pressure was already set once. Use the .update_pressure(self) method to update pressure based on current level.') def set_influx(self,influx): # sets influx to the reservoir in m³/s @@ -143,7 +149,7 @@ class Ausgleichsbecken_class: ss_outflux = ss_influx ss_influx_vel = abs(ss_influx/self.area) ss_outflux_vel = abs(ss_outflux/self.area_out) - # see confluence doc for explaination on how to arrive at the ss pressure formula + # see word document 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) # use setter methods to set the attributes to their steady state values @@ -155,6 +161,7 @@ class Ausgleichsbecken_class: # getter - return attributes def get_info(self, full = False): # prints out the info on the current state of the reservoir + # full = True gives more info new_line = '\n' if self.pressure != np.inf: p = pressure_conversion(self.pressure,self.pressure_unit,self.pressure_unit_disp) @@ -189,7 +196,8 @@ class Ausgleichsbecken_class: 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"----------------------------- {new_line}") - + + # print the info to console print(print_str) def get_current_influx(self): @@ -210,12 +218,15 @@ class Ausgleichsbecken_class: # update methods - update attributes based on some parameter def update_level(self,timestep,set_flag=False): # 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 then convert the new volume to a level by assuming a cuboid reservoir + # there is no call of the update_volume() function because I need the updated level from half a timestep in the reservoir evolution + # if update_volume() was called within this function, the script would produce wrong results. net_flux = self.influx-self.outflux delta_level = net_flux*timestep/self.area level_new = (self.level+delta_level) - if level_new < 0.1: - raise Exception('Ausgleichsbecken leer') + # raise exception error if level in reservoir falls below 0.01 ######################### has to be commented out if used in loop + if level_new < 0.01: + raise Exception('Reservoir ran emtpy') # set flag is necessary because update_level() is used to get a halfstep value in the time evoultion if set_flag == True: self.set_level(level_new,display_warning=False) @@ -224,7 +235,7 @@ class Ausgleichsbecken_class: def update_pressure(self,set_flag=False): # update pressure based on level and flux velocities - # see confluence doc for explaination + # see word document for explaination influx_vel = abs(self.influx/self.area) outflux_vel = abs(self.outflux/self.area_out) p_new = self.density*self.g*self.level+self.density*outflux_vel*(influx_vel-outflux_vel) @@ -245,6 +256,7 @@ class Ausgleichsbecken_class: #methods def timestep_reservoir_evolution(self): # update outflux, level, pressure and volume based on current pipeline pressure and waterlevel in reservoir + # solve the FODE of the outflux velocity for one timestep using explicit four step Runge-Kutta method # get some variables dt = self.timestep diff --git a/Druckrohrleitung/Druckrohrleitung_class_file.py b/Druckrohrleitung/Druckrohrleitung_class_file.py index 954559b..9175e18 100644 --- a/Druckrohrleitung/Druckrohrleitung_class_file.py +++ b/Druckrohrleitung/Druckrohrleitung_class_file.py @@ -1,5 +1,8 @@ -import os -import sys +# import modules for general use +import os # to import functions from other folders +import sys # to import functions from other folders +from logging import \ + exception # to throw an exception when a specific condition is met import numpy as np @@ -39,6 +42,7 @@ class Druckrohrleitung_class: g = 9.81 # m/s² gravitational acceleration # init + # see docstring below def __init__(self,total_length,diameter,pipeline_head,number_segments,Darcy_friction_factor,pw_vel,timestep,pressure_unit_disp,rho=1000): """ Creates a reservoir with given attributes in this order: \n @@ -154,6 +158,7 @@ class Druckrohrleitung_class: 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 + # dynamic pressure does not play a role, because it has the same influence on both sides of the equation (constant flow velocity) and therefore cancels out 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) # set the initial conditions @@ -162,6 +167,7 @@ class Druckrohrleitung_class: # getter - return attributes def get_info(self): + # prints out the info on the current state of the reservoir new_line = '\n' angle_deg = round(self.angle/np.pi*180,3) @@ -182,8 +188,11 @@ class Druckrohrleitung_class: 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"----------------------------- {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 via the {new_line} \ + get_current_velocity_distribution() and get_current_pressure_distribution() methods of the pipeline object. {new_line} \ + See also get_lowest_XXX_per_node() and get_highest_XXX_per_node() methods.") + # print the info to console print(print_str) def get_current_pressure_distribution(self,disp_flag=False): @@ -200,12 +209,14 @@ class Druckrohrleitung_class: return self.v*self.A def get_lowest_pressure_per_node(self,disp_flag=False): + # 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_min,self.pressure_unit,self.pressure_unit_disp) elif disp_flag == False: # stay in Pa return self.p_min def get_highest_pressure_per_node(self,disp_flag=False): + # 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_max,self.pressure_unit,self.pressure_unit_disp) elif disp_flag == False: # stay in Pa @@ -244,7 +255,7 @@ class Druckrohrleitung_class: g = self.g # graviational acceleration alpha = self.angle # pipeline angle - # Vectorize this loop? + # Vectorized loop see below for i in range(1,nn-1): self.v[i] = 0.5*(self.v_old[i+1]+self.v_old[i-1])-0.5/(rho*c)*(self.p_old[i+1]-self.p_old[i-1]) \ +dt*g*np.sin(alpha)-f_D*dt/(4*D)*(abs(self.v_old[i+1])*self.v_old[i+1]+abs(self.v_old[i-1])*self.v_old[i-1]) @@ -265,6 +276,7 @@ class Druckrohrleitung_class: self.v_old = self.v.copy() def timestep_characteristic_method_vectorized(self): + # faster then above # 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 diff --git a/Kraftwerk/Kraftwerk_class_file.py b/Kraftwerk/Kraftwerk_class_file.py index aa1ad86..b0a57e1 100644 --- a/Kraftwerk/Kraftwerk_class_file.py +++ b/Kraftwerk/Kraftwerk_class_file.py @@ -1,20 +1,33 @@ +# import modules for general use +import os # to import functions from other folders +import sys # to import functions from other folders +from logging import \ + exception # to throw an exception when a specific condition is met + import numpy as np -#importing Druckrohrleitung -import sys -import os -current = os.path.dirname(os.path.realpath('Main_Programm.ipynb')) + +#importing pressure conversion function +current = os.path.dirname(os.path.realpath(__file__)) parent = os.path.dirname(current) sys.path.append(parent) from functions.pressure_conversion import pressure_conversion from Turbinen.Turbinen_class_file import Francis_Turbine + class Kraftwerk_class: g = 9.81 def __init__(self): + # create an empty powerhouse + # see add_turbine() method self.turbines = [] self.n_turbines = 0 + def add_turbine(self,turbine): + # add a turbine object from the turbine class + self.turbines.append(turbine) + self.n_turbines += 1 + # setter def set_LAs(self,LA_vec,display_warning=True): for i in range(self.n_turbines): @@ -61,14 +74,11 @@ class Kraftwerk_class: # methods def identify_Q_proportion(self): + # calculate the proportions of the nominal fluxes of all turbines in the powerhouse Q_n_vec = np.zeros(self.n_turbines) for i in range(self.n_turbines): Q_n_vec[i] = self.turbines[i].get_Q_n() self.Q_prop = Q_n_vec/np.sum(Q_n_vec) - - def add_turbine(self,turbine): - self.turbines.append(turbine) - self.n_turbines += 1 def update_LAs(self,LA_soll_vec): for i in range(self.n_turbines): @@ -77,7 +87,7 @@ class Kraftwerk_class: 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 # new pressure from the forward characteristic are not perfectly compatible. - # Therefore, iterate the flux and the pressure so long, until they converge + # Therefore, iterate the flux and the pressure so long, until they converge - i honestly have no idea why that works :D (steady state test prove it right ¯\_(ツ)_/¯) eps = 1e-12 # convergence criterion: iteration change < eps iteration_change = 1. # change in Q from one iteration to the next diff --git a/Regler/Regler_class_file.py b/Regler/Regler_class_file.py index 5b9135b..01e389a 100644 --- a/Regler/Regler_class_file.py +++ b/Regler/Regler_class_file.py @@ -1,4 +1,5 @@ import numpy as np + #based on https://en.wikipedia.org/wiki/PID_controller#Discrete_implementation # performance parameters for controllers @@ -129,7 +130,7 @@ class PI_controller_class: def get_info(self): new_line = '\n' # :<10 pads the self.value to be 10 characters wide - print_str = (f"Turbine has the following attributes: {new_line}" + print_str = (f"Controller has the following attributes: {new_line}" f"----------------------------- {new_line}" f"Type = PI Controller {new_line}" f"Setpoint = {self.SP:<10} {new_line}" diff --git a/Turbinen/Turbinen_class_file.py b/Turbinen/Turbinen_class_file.py index b39a01e..2958db2 100644 --- a/Turbinen/Turbinen_class_file.py +++ b/Turbinen/Turbinen_class_file.py @@ -163,7 +163,7 @@ class Francis_Turbine: 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 # new pressure from the forward characteristic are not perfectly compatible. - # Therefore, iterate the flux and the pressure so long, until they converge + # Therefore, iterate the flux and the pressure so long, until they converge - i honestly have no idea why that works :D (steady state test prove it right ¯\_(ツ)_/¯) eps = 1e-12 # convergence criterion: iteration change < eps iteration_change = 1. # change in Q from one iteration to the next