diff --git a/Ausgleichsbecken/Ausgleichsbecken_class_file.py b/Ausgleichsbecken/Ausgleichsbecken_class_file.py index b00d957..0092846 100644 --- a/Ausgleichsbecken/Ausgleichsbecken_class_file.py +++ b/Ausgleichsbecken/Ausgleichsbecken_class_file.py @@ -1,3 +1,4 @@ +from logging import exception import numpy as np #importing pressure conversion function @@ -8,8 +9,19 @@ parent = os.path.dirname(current) sys.path.append(parent) from functions.pressure_conversion import pressure_conversion -def FODE_function(x, h, alpha, p, rho=1000., g=9.81): - f = x*abs(x)/h*alpha+g-p/(rho*h) +def FODE_function(x,h,A,A_a,p,rho,g): + # (FODE ... first order differential equation) + # 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 + # x ... effusion velocity + # h ... level in the reservoir + # A_a ... Outflux_Area + # A ... Reservoir_Area + # g ... gravitational acceleration + # rho ... density of the liquid in the reservoir + f = x*abs(x)/h*(A_a/A-1)+g-p/(rho*h) return f @@ -19,67 +31,84 @@ class Ausgleichsbecken_class: # units are used to label graphs and print units are used to have a bearable format 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' time_unit = 's' velocity_unit = r'$\mathrm{m}/\mathrm{s}$' volume_unit = r'$\mathrm{m}^3$' area_unit_print = 'm²' area_outflux_unit_print = 'm²' + density_unit_print = 'kg/m³' flux_unit_print = 'm³/s' level_unit_print = 'm' time_unit_print = 's' + pressure_unit_print = '--' # will be set by .set_pressure() method velocity_unit_print = 'm/s' volume_unit_print = 'm³' - g = 9.81 # m/s² - rho = 1000 # kg/m³ + g = 9.81 # m/s² gravitational acceleration + # init - def __init__(self,area,outflux_area,level_min = 0,level_max = np.inf ,timestep = 1): + def __init__(self,area,outflux_area,level_min = 0,level_max = np.inf ,timestep = 1,rho = 1000): self.area = area # base area of the rectangular structure self.area_outflux = outflux_area # area of the outlet towards the pipeline + self.density = rho # density of the liquid in the system self.level_min = level_min # lowest allowed water level self.level_max = level_max # highest allowed water level self.timestep = timestep # timestep of the simulation # initialize for get_info - self.level = "--" self.influx = "--" + self.level = "--" self.outflux = "--" self.volume = "--" # setter - def set_volume(self): + def update_volume(self): + # sets volume in reservoir based on self.level self.volume = self.level*self.area def set_initial_level(self,initial_level): - self.level = initial_level - self.set_volume() + # sets the level in the reservoir and should only be called during initialization + if self.level == '--': + self.level = initial_level + self.update_volume() + else: + raise Exception('Initial level was already set once. Use the .update_level(self,timestep) method to update level based on net flux.') def set_influx(self,influx): + # sets influx to the reservoir in m³/s + # positive influx means that liquid flows into the reservoir self.influx = influx def set_outflux(self,outflux): - self.outflux = outflux - self.outflux_vel = outflux/self.area_outflux + # sets outflux to the reservoir in m³/s + # positive outflux means that liquid flows out of reservoir the reservoir + self.outflux = outflux - def set_pressure(self,pressure,pressure_unit,display_pressure_unit): + def set_pressure(self,pressure,display_pressure_unit): + # sets the static pressure present at the outlet of the reservoir + # units are used to convert and display the pressure self.pressure = pressure - self.pressure_unit = pressure_unit self.pressure_unit_print = display_pressure_unit - def set_steady_state(self,ss_influx,ss_level,pressure_unit,display_pressure_unit): + def set_steady_state(self,ss_influx,ss_level,display_pressure_unit): + # find the steady state (ss) condition in which the net flux is zero + # set pressure acting on the outflux so that the level stays constant ss_outflux = ss_influx ss_outflux_vel = ss_outflux/self.area_outflux - ss_pressure = self.rho*self.g*ss_level-ss_outflux_vel**2*self.rho/2 + ss_pressure = self.density*self.g*ss_level-ss_outflux_vel**2*self.density/2 self.set_initial_level(ss_level) self.set_influx(ss_influx) self.set_outflux(ss_outflux) - self.set_pressure(ss_pressure,pressure_unit,display_pressure_unit) + self.set_pressure(ss_pressure,display_pressure_unit) + # getter def get_info(self, full = False): new_line = '\n' @@ -101,6 +130,7 @@ class Ausgleichsbecken_class: f"Current outflux vel = {round(self.outflux_vel,3):<10} {self.velocity_unit_print} {new_line}" f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}" f"Simulation timestep = {self.timestep:<10} {self.time_unit_print} {new_line}" + f"Density of liquid = {self.density:<10} {self.density_unit_print} {new_line}" f"----------------------------- {new_line}") else: # :<10 pads the self.value to be 10 characters wide @@ -116,9 +146,27 @@ class Ausgleichsbecken_class: print(print_str) + def get_current_level(self): + return self.level + + def get_current_influx(self): + return self.influx + + def get_current_outflux(self): + return self.outflux + + def get_current_volume(self): + return self.volume + + def get_current_pressure(self): + return self.pressure + + # methods def update_level(self,timestep): + # 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 net_flux = self.influx-self.outflux delta_V = net_flux*timestep new_level = (self.volume+delta_V)/self.area @@ -127,21 +175,25 @@ class Ausgleichsbecken_class: def e_RK_4(self): # update outflux and outflux velocity based on current pipeline pressure and waterlevel in reservoir - yn = self.outflux_vel + yn = self.outflux/self.area_outflux # outflux velocity h = self.level dt = self.timestep p = self.pressure - # assume constant pipeline pressure during timestep (see comments in main_programm) + # assume constant pipeline pressure during timestep + # e_RK_4 timestep is way smalle than timestep of characteristic method, so this should be a valid approx. + # (furthermore I have no idea how to approximate p_hs otherwise :/ ) p_hs = self.pressure - alpha = (self.area_outflux/self.area-1) + A_a = self.area_outflux + A = self.area h_hs = self.update_level(dt/2) + rho = self.density + g = self.g # explicit 4 step Runge Kutta Y1 = yn - Y2 = yn + dt/2*FODE_function(Y1, h, alpha, self.pressure) - Y3 = yn + dt/2*FODE_function(Y2, h_hs, alpha, p_hs) - Y4 = yn + dt*FODE_function(Y3, h_hs, alpha, p_hs) - ynp1 = yn + dt/6*(FODE_function(Y1, h, alpha, p)+2*FODE_function(Y2, h_hs, alpha, p_hs)+ \ - 2*FODE_function(Y3, h_hs, alpha, p_hs)+ FODE_function(Y4, h, alpha, p)) + Y2 = yn + dt/2*FODE_function(Y1,h,A,A_a,self.pressure,rho,g) + Y3 = yn + dt/2*FODE_function(Y2,h_hs,A,A_a,p_hs,rho,g) + Y4 = yn + dt*FODE_function(Y3,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)) - self.outflux_vel = ynp1 self.outflux = ynp1*self.area_outflux diff --git a/Ausgleichsbecken/Ausgleichsbecken_test.ipynb b/Ausgleichsbecken/Ausgleichsbecken_test.ipynb index eea80ce..4f72ad4 100644 --- a/Ausgleichsbecken/Ausgleichsbecken_test.ipynb +++ b/Ausgleichsbecken/Ausgleichsbecken_test.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 12, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -21,7 +21,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -46,7 +46,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -76,7 +76,7 @@ " V.pressure = pressure_vec[i]\n", " V.e_RK_4()\n", " V.level = V.update_level(V.timestep)\n", - " V.set_volume()\n", + " V.update_volume()\n", " outflux_vec[i+1] = V.outflux\n", " level_vec[i+1] = V.level\n", " if V.level < total_min_level:\n", @@ -87,7 +87,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -141,7 +141,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3.8.13 ('DT_Slot_3')", + "display_name": "Python 3.8.13 ('Georg_DT_Slot3')", "language": "python", "name": "python3" }, @@ -160,7 +160,7 @@ "orig_nbformat": 4, "vscode": { "interpreter": { - "hash": "4a28055eb8a3160fa4c7e4fca69770c4e0a1add985300856aa3fcf4ce32a2c48" + "hash": "84fb123bdc47ab647d3782661abcbe80fbb79236dd2f8adf4cef30e8755eb2cd" } } }, diff --git a/functions/pressure_conversion.py b/functions/pressure_conversion.py index 96744c4..d4a91f9 100644 --- a/functions/pressure_conversion.py +++ b/functions/pressure_conversion.py @@ -32,7 +32,7 @@ def pa_to_atm(p): def pa_to_psi(p): return p/6894.8 -def pressure_conversion(pressure, input_unit = 'bar', target_unit = 'Pa'): +def pressure_conversion(pressure, input_unit = 'bar', target_unit = 'Pa', return_unit = False): p = pressure if input_unit.lower() == 'bar': p_pa = bar_to_pa(p) @@ -50,20 +50,27 @@ def pressure_conversion(pressure, input_unit = 'bar', target_unit = 'Pa'): raise Exception('Given input unit not recognised. \n Known units are: Pa, bar, mWs, Torr, atm, psi') if target_unit.lower() == 'bar': - return pa_to_bar(p_pa), target_unit + return_vec = [pa_to_bar(p_pa), target_unit] elif target_unit.lower() == 'mws': - return pa_to_mWS(p_pa), target_unit + return_vec = [pa_to_mWS(p_pa), target_unit] elif target_unit.lower() == 'torr': - return pa_to_torr(p_pa), target_unit + return_vec = [pa_to_torr(p_pa), target_unit] elif target_unit.lower() == 'atm': - return pa_to_atm(p_pa), target_unit + return_vec = [pa_to_atm(p_pa), target_unit] elif target_unit.lower() =='psi': - return pa_to_psi(p_pa), target_unit + return_vec = [pa_to_psi(p_pa), target_unit] elif target_unit.lower() == 'pa': - return p_pa, target_unit + return_vec = [p_pa, target_unit] else: raise Exception('Given target unit not recognised. \n Known units are: Pa, bar, mWs, Torr, atm, psi') + if return_unit == True: + # return with pressure unit + return return_vec + else: + # return without pressure unit + return return_vec[0] + # testing_pressure_conversion if __name__ == '__main__': p = 1 @@ -72,6 +79,6 @@ if __name__ == '__main__': for input_unit in unit_dict: for target_unit in unit_dict: - converted_p = pressure_conversion(p,input_unit,target_unit) + converted_p = pressure_conversion(p,input_unit,target_unit,return_unit=False) print(input_unit,target_unit) print(converted_p) \ No newline at end of file