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

added comments in preparation of merge

Browse Source
This commit is contained in:
Brantegger Georg
2023-02-06 09:57:23 +01:00
parent 1b70fbdab8
commit f855799ec1
5 changed files with 67 additions and 32 deletions
Download Patch File Download Diff File Expand all files Collapse all files

42
Ausgleichsbecken/Ausgleichsbecken_class_file.py
View File

@@ -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
# 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)
@@ -190,6 +197,7 @@ class Ausgleichsbecken_class:
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

20
Druckrohrleitung/Druckrohrleitung_class_file.py
View File

@@ -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

28
Kraftwerk/Kraftwerk_class_file.py
View File

@@ -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,15 +74,12 @@ 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):
self.turbines[i].update_LA(LA_soll_vec[i])
@@ -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

3
Regler/Regler_class_file.py
View File

@@ -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}"

2
Turbinen/Turbinen_class_file.py
View File

@@ -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