From ba696444bb0da62d5e0170595c7261fabea9e2b4 Mon Sep 17 00:00:00 2001 From: Brantegger Georg Date: Wed, 3 Aug 2022 15:56:56 +0200 Subject: [PATCH] fix for numerical runaway of rounding errors due to turbine-pipeline interatction via a convergence method in the turbine and a "damping" trick on the reservoir velocity plus: code cleanup with consistent naming of variables --- .../Ausgleichsbecken_class_file.py | 217 +++++---- .../Ausgleichsbecken_test_steady_state.ipynb | 207 +++++---- .../Druckrohrleitung_class_file.py | 143 +++--- .../Druckrohrleitung_test_steady_state.ipynb | 176 ++++---- Regler/Pegelregler_test.ipynb | 285 +++++------- Regler/Regler_class_file.py | 6 +- Turbinen/Turbinen_class_file.py | 131 ++++-- Turbinen/Turbinen_test_steady_state.ipynb | 370 ++++++++++++++++ Turbinen/{ => old}/convergence_turbine.py | 3 +- .../{ => old}/turbine_convergence_test.ipynb | 120 +++-- Untertweng.ipynb | 360 --------------- Untertweng_mit_Pegelregler.ipynb | 415 ++++++++++-------- untertweng.txt | 22 - 13 files changed, 1257 insertions(+), 1198 deletions(-) create mode 100644 Turbinen/Turbinen_test_steady_state.ipynb rename Turbinen/{ => old}/convergence_turbine.py (99%) rename Turbinen/{ => old}/turbine_convergence_test.ipynb (75%) delete mode 100644 Untertweng.ipynb delete mode 100644 untertweng.txt diff --git a/Ausgleichsbecken/Ausgleichsbecken_class_file.py b/Ausgleichsbecken/Ausgleichsbecken_class_file.py index 89faa97..2eab201 100644 --- a/Ausgleichsbecken/Ausgleichsbecken_class_file.py +++ b/Ausgleichsbecken/Ausgleichsbecken_class_file.py @@ -15,10 +15,10 @@ def FODE_function(x_out,h,A,A_a,p,rho,g): # https://www.youtube.com/watch?v=8HO2LwqOhqQ # adapted for a pressurized pipeline into which the reservoir effuses # and flow direction - # x_out ... effusion velocity + # x_out ... effusion velocity # h ... level in the reservoir - # A_a ... Outflux_Area - # A ... Reservoir_Area + # A_a ... Area_outflux + # A ... Area_reservoir_base # g ... gravitational acceleration # rho ... density of the liquid in the reservoir f = x_out*abs(x_out)/h*(A_a/A-1.)+g-p/(rho*h) @@ -28,168 +28,202 @@ def FODE_function(x_out,h,A,A_a,p,rho,g): class Ausgleichsbecken_class: # units # make sure that units and print units are the same - # units are used to label graphs and print units are used to have a bearable format when using pythons print() - 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$' + # 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_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' - pressure_unit_print = '--' # will be set by .set_pressure() method - time_unit_print = 's' - velocity_unit_print = 'm/s' - volume_unit_print = 'm³' + area_unit_disp = 'm²' + area_outflux_unit_disp = 'm²' + density_unit_disp = 'kg/m³' + flux_unit_disp = 'm³/s' + level_unit_disp = 'm' + time_unit_disp = 's' + velocity_unit_disp = 'm/s' + volume_unit_disp = 'm³' g = 9.81 # m/s² gravitational acceleration # init - 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 + def __init__(self,area,area_outflux,timestep,pressure_unit_disp,level_min=0,level_max=np.inf,rho = 1000.): + self.area = area # base area of the cuboid reservoir + self.area_out = area_outflux # 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.pressure_unit_disp = pressure_unit_disp # pressure unit for displaying + self.timestep = timestep # timestep in the time evolution method # initialize for get_info - self.influx = "--" - self.level = "--" - self.outflux = "--" - self.volume = "--" + self.influx = "--" + self.outflux = "--" + self.level = "--" + self.pressure = "--" + self.volume = "--" # setter def set_initial_level(self,initial_level): - # sets the level in the reservoir and should only be called during initialization + # sets the initial level in the reservoir and should only be called during initialization if self.level == '--': self.level = initial_level + self.update_volume(set_flag=True) else: - raise Exception('Initial level was already set once. Use the .update_level(self,timestep) method to update level based on net flux.') + raise Exception('Initial level was already set once. Use the .update_level(self,timestep,set_flag=True) method to update level based on net flux.') - def set_level(self,level): - self.level = level + 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 + if self.pressure == '--': + self.pressure = initial_pressure + else: + raise Exception('Initial pressure was already set once. Use the .update_pressure(self) method to update pressure based current level.') def set_influx(self,influx): # 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): + def set_outflux(self,outflux,display_warning=True): # sets outflux to the reservoir in m³/s # positive outflux means that liquid flows out of reservoir the reservoir + if display_warning == True: + print('You are setting the outflux from the reservoir manually. \n \ + This is not an intended use of this method. \n \ + Refer to the timestep_reservoir_evolution() method instead.') self.outflux = outflux - def set_initial_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_print = display_pressure_unit + def set_level(self,level,display_warning=True): + # sets level in the reservoir in m + if display_warning == True: + print('You are setting the level of the reservoir manually. \n \ + This is not an intended use of this method. \n \ + Refer to the update_level() method instead.') + self.level = level - def set_pressure(self,pressure): - # sets the static pressure present at the outlet of the reservoir - self.pressure = pressure + def set_pressure(self,pressure,display_warning=True): + # sets pressure in the pipeline just below the reservoir in Pa + if display_warning == True: + print('You are setting the pressure below the reservoir manually. \n \ + This is not an intended use of this method. \n \ + Refer to the update_pressure() method instead.') + self.pressure = pressure - def set_steady_state(self,ss_influx,ss_level,display_pressure_unit): + def set_volume(self,volume,display_warning=True): + if display_warning == True: + print('You are setting the volume in the reservoir manually. \n \ + This is not an intended use of this method. \n \ + Refer to the .update_volume() or set_initial_level() method instead.') + self.volume = volume + + def set_steady_state(self,ss_influx,ss_level): # set the steady state (ss) condition in which the net flux is zero # set pressure acting on the outflux area so that the level stays constant ss_outflux = ss_influx ss_influx_vel = abs(ss_influx/self.area) - ss_outflux_vel = abs(ss_outflux/self.area_outflux) + ss_outflux_vel = abs(ss_outflux/self.area_out) ss_pressure = self.density*self.g*ss_level+self.density*ss_outflux_vel*(ss_influx_vel-ss_outflux_vel) self.set_influx(ss_influx) self.set_initial_level(ss_level) - self.set_initial_pressure(ss_pressure,display_pressure_unit) - self.set_outflux(ss_outflux) + self.set_initial_pressure(ss_pressure) + self.set_outflux(ss_outflux,display_warning=False) # getter def get_info(self, full = False): new_line = '\n' - p = pressure_conversion(self.pressure,self.pressure_unit,self.pressure_unit_print) - outflux_vel = self.outflux/self.area_outflux + p = pressure_conversion(self.pressure,self.pressure_unit,self.pressure_unit_disp) + outflux_vel = self.outflux/self.area_out if full == True: # :<10 pads the self.value to be 10 characters wide print_str = (f"The cuboid reservoir has the following attributes: {new_line}" f"----------------------------- {new_line}" - f"Base area = {self.area:<10} {self.area_unit_print} {new_line}" - f"Outflux area = {round(self.area_outflux,3):<10} {self.area_outflux_unit_print} {new_line}" - f"Current level = {self.level:<10} {self.level_unit_print}{new_line}" - f"Critical level low = {self.level_min:<10} {self.level_unit_print} {new_line}" - f"Critical level high = {self.level_max:<10} {self.level_unit_print} {new_line}" - f"Volume in reservoir = {self.volume:<10} {self.volume_unit_print} {new_line}" - f"Current influx = {self.influx:<10} {self.flux_unit_print} {new_line}" - f"Current outflux = {self.outflux:<10} {self.flux_unit_print} {new_line}" - f"Current outflux vel = {round(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"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"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 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"Current influx = {self.influx:<10} {self.flux_unit_disp} {new_line}" + f"Current outflux = {self.outflux:<10} {self.flux_unit_disp} {new_line}" + f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_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"Density of liquid = {self.density:<10} {self.density_unit_disp} {new_line}" f"----------------------------- {new_line}") else: # :<10 pads the self.value to be 10 characters wide print_str = (f"The current attributes are: {new_line}" f"----------------------------- {new_line}" - f"Current level = {self.level:<10} {self.level_unit_print}{new_line}" - f"Volume in reservoir = {self.volume:<10} {self.volume_unit_print} {new_line}" - f"Current influx = {self.influx:<10} {self.flux_unit_print} {new_line}" - f"Current outflux = {self.outflux:<10} {self.flux_unit_print} {new_line}" - f"Current outflux vel = {round(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"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 influx = {self.influx:<10} {self.flux_unit_disp} {new_line}" + f"Current outflux = {self.outflux:<10} {self.flux_unit_disp} {new_line}" + f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_disp} {new_line}" + f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_disp} {new_line}" f"----------------------------- {new_line}") 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_level(self): + return self.level def get_current_pressure(self): return self.pressure + def get_current_volume(self): + return self.volume - -# methods - - def update_level(self,timestep): +# update methods + 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 - - # cannot set new level directly in this method, because it gets called to calcuate during the Runge Kutta - # to calculate a ficticious level at half the timestep net_flux = self.influx-self.outflux delta_level = net_flux*timestep/self.area - new_level = (self.level+delta_level) - return new_level + level_new = (self.level+delta_level) + if set_flag == True: + self.set_level(level_new,display_warning=False) + elif set_flag == False: + return level_new - def update_pressure(self): + def update_pressure(self,set_flag=False): influx_vel = abs(self.influx/self.area) - outflux_vel = abs(self.outflux/self.area_outflux) + outflux_vel = abs(self.outflux/self.area_out) p_new = self.density*self.g*self.level+self.density*outflux_vel*(influx_vel-outflux_vel) - return p_new + if set_flag ==True: + self.set_pressure(p_new,display_warning=False) + elif set_flag == False: + return p_new + def update_volume(self,set_flag=False): + volume_new = self.level*self.area + if set_flag == True: + self.set_volume(volume_new,display_warning=False) + elif set_flag == False: + return volume_new + +#methods def timestep_reservoir_evolution(self): # update outflux and outflux velocity based on current pipeline pressure and waterlevel in reservoir dt = self.timestep rho = self.density g = self.g A = self.area - A_a = self.area_outflux + A_a = self.area_out yn = self.outflux/A_a # outflux velocity h = self.level h_hs = self.update_level(dt/2) @@ -203,10 +237,7 @@ class Ausgleichsbecken_class: ynp1 = yn + dt/6*(FODE_function(Y1,h,A,A_a,p,rho,g)+2*FODE_function(Y2,h_hs,A,A_a,p_hs,rho,g)+ \ 2*FODE_function(Y3,h_hs,A,A_a,p_hs,rho,g)+ FODE_function(Y4,h,A,A_a,p,rho,g)) - new_outflux = ynp1*A_a - new_level = self.update_level(dt) - new_pressure = self.update_pressure() - - self.set_outflux(new_outflux) - self.set_level(new_level) - self.set_pressure(new_pressure) \ No newline at end of file + self.set_outflux(ynp1*A_a,display_warning=False) + self.update_level(dt,set_flag=True) + self.update_volume(set_flag=True) + self.update_pressure(set_flag=True) diff --git a/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb b/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb index e420f63..918aa20 100644 --- a/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb +++ b/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 7, + "execution_count": 29, "metadata": {}, "outputs": [], "source": [ @@ -21,151 +21,142 @@ }, { "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [], - "source": [ - "L = 1000.\n", - "n = 10000 # number of pipe segments in discretization\n", - "c = 400. \n", - "dx = L/n # length of each pipe segment\n", - "dt = dx/c \n", - "\n", - "# # define constants\n", - "# initial_level = 10.1 # m\n", - "# initial_influx = 0.8 # m³/s\n", - "# conversion_pressure_unit = 'mWS'\n", - "\n", - "# area_base = 75. # m²\n", - "# area_outflux = (0.9/2)**2*np.pi # m²\n", - "# critical_level_low = 0. # m\n", - "# critical_level_high = 10. # m\n", - "# simulation_timestep = dt # s\n", - "\n", - "# # for while loop\n", - "# total_min_level = 0.01 # m\n", - "# total_max_time = 100 # s\n", - "\n", - "# nt = int(total_max_time//simulation_timestep)" - ] - }, - { - "cell_type": "code", - "execution_count": 9, + "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "# define constants\n", - "initial_level = 10.1 # m\n", - "initial_influx = 1. # m³/s\n", - "# initial_outflux = 1. # m³/s\n", - "# initial_pipeline_pressure = 10.\n", - "# initial_pressure_unit = 'mWS'\n", - "conversion_pressure_unit = 'mWS'\n", "\n", - "area_base = 75. # m²\n", - "area_outflux = 2. # m²\n", - "critical_level_low = 0. # m\n", - "critical_level_high = 10. # m\n", - "simulation_timestep = dt # s\n", + " # for physics\n", + "g = 9.81 # [m/s²] gravitational acceleration \n", + "rho = 1000. # [kg/m³] density of water \n", + "pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", + "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n", "\n", - "# for while loop\n", - "total_min_level = 0.01 # m\n", - "total_max_time = 1000 # s\n", "\n", - "nt = int(total_max_time//simulation_timestep)" + " # for Turbine\n", + "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", + "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", + "Tur_closingTime = 90. # [s] closing time of turbine\n", + "\n", + "\n", + " # for PI controller\n", + "Con_targetLevel = 8. # [m]\n", + "Con_K_p = 0.1 # [-] proportional constant of PI controller\n", + "Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\n", + "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n", + "\n", + "\n", + " # for pipeline\n", + "Pip_length = (535.+478.) # [m] length of pipeline\n", + "Pip_dia = 0.9 # [m] diameter of pipeline\n", + "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", + "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", + "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", + "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", + "Pip_f_D = 0.014 # [-] Darcy friction factor\n", + "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", + " # derivatives of the pipeline constants\n", + "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", + "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", + "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", + "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", + "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", + "\n", + "\n", + " # for reservoir\n", + "Res_area_base = 5. # [m²] total base are of the cuboid reservoir \n", + "Res_area_out = Pip_area # [m²] outflux area of the reservoir, given by pipeline area\n", + "Res_level_crit_lo = 0. # [m] for yet-to-be-implemented warnings\n", + "Res_level_crit_hi = np.inf # [m] for yet-to-be-implemented warnings\n", + "Res_dt_approx = 1e-3 # [s] approx. timestep of reservoir time evolution to ensure numerical stability (see Res_nt why approx.)\n", + "Res_nt = max(1,int(Pip_dt//Res_dt_approx)) # [1] number of timesteps of the reservoir time evolution within one timestep of the pipeline\n", + "Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n", + "\n", + " # for general simulation\n", + "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", + "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", + "simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", + "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", + "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 31, "metadata": {}, "outputs": [], "source": [ - "%matplotlib qt\n", + "# create objects\n", "\n", - "V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,simulation_timestep)\n", - "# V.set_initial_level(initial_level) \n", - "# V.set_influx(initial_influx)\n", - "# V.set_outflux(initial_outflux)\n", - "# V.set_initial_pressure(pressure_conversion(initial_pipeline_pressure,input_unit = initial_pressure_unit, target_unit = 'Pa'),conversion_pressure_unit)\n", - "# V.pressure = converted_pressure\n", - "V.set_steady_state(initial_influx,initial_level,conversion_pressure_unit)\n", + "# Upstream reservoir\n", + "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n", + "reservoir.set_steady_state(flux_init,level_init)\n", "\n", - "time_vec = np.arange(0,nt+1,1)*simulation_timestep\n", - "outflux_vec = np.zeros_like(time_vec)\n", - "outflux_vec[0] = V.get_current_outflux()\n", - "level_vec = np.zeros_like(time_vec)\n", - "level_vec[0] = V.get_current_level()\n", - "pressure_vec = np.zeros_like(time_vec)\n", - "pressure_vec[0] = V.get_current_pressure()\n", - "\n", - "# pressure_vec = np.full_like(time_vec,converted_pressure)*((np.sin(time_vec)+1)*np.exp(-time_vec/50))\n", - " \n", - "i_max = -1\n", + "reservoir.get_info(full=True)\n", "\n", + "# initialize vectors\n", + "outflux_vec = np.zeros_like(t_vec)\n", + "outflux_vec[0] = reservoir.get_current_outflux()\n", + "level_vec = np.zeros_like(t_vec)\n", + "level_vec[0] = reservoir.get_current_level()\n", + "volume_vec = np.zeros_like(t_vec)\n", + "volume_vec[0] = reservoir.get_current_volume()\n", + "pressure_vec = np.zeros_like(t_vec)\n", + "pressure_vec[0] = reservoir.get_current_pressure()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# time loop\n", "for i in range(1,nt+1):\n", - " V.set_pressure(pressure_vec[i-1])\n", - " V.set_outflux(outflux_vec[i-1])\n", - " V.timestep_reservoir_evolution()\n", - " outflux_vec[i] = V.get_current_outflux()\n", - " level_vec[i] = V.get_current_level()\n", - " pressure_vec[i] = V.get_current_pressure()\n", - " if V.level < total_min_level:\n", - " i_max = i\n", - " break\n", - "\n" + " # if i == 500:\n", + " # reservoir.set_influx(0.)\n", + " reservoir.set_pressure(pressure_vec[i-1],display_warning=False)\n", + " reservoir.set_outflux(outflux_vec[i-1],display_warning=False)\n", + " for it_res in range(Res_nt):\n", + " reservoir.timestep_reservoir_evolution() \n", + " \n", + " outflux_vec[i] = reservoir.get_current_outflux()\n", + " level_vec[i] = reservoir.get_current_level()\n", + " pressure_vec[i] = reservoir.get_current_pressure()\n", + "\n", + " reservoir.get_info()" ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 32, "metadata": {}, "outputs": [], "source": [ - "\n", + "%matplotlib qt5\n", "fig1, (ax1, ax2, ax3) = plt.subplots(3, 1)\n", "fig1.set_figheight(10)\n", "fig1.suptitle('Ausgleichsbecken')\n", "\n", - "ax1.plot(time_vec[:i_max],level_vec[:i_max], label='Water level')\n", - "ax1.set_ylabel(r'$h$ ['+V.level_unit+']')\n", - "ax1.set_xlabel(r'$t$ ['+V.time_unit+']')\n", + "ax1.plot(t_vec,level_vec, label='Water level')\n", + "ax1.set_ylabel(r'$h$ ['+reservoir.level_unit+']')\n", + "ax1.set_xlabel(r'$t$ ['+reservoir.time_unit+']')\n", "ax1.legend()\n", "\n", - "ax2.plot(time_vec[:i_max],outflux_vec[:i_max], label='Outflux')\n", - "ax2.set_ylabel(r'$Q_{out}$ ['+V.flux_unit+']')\n", - "ax2.set_xlabel(r'$t$ ['+V.time_unit+']')\n", + "ax2.plot(t_vec,outflux_vec, label='Outflux')\n", + "ax2.set_ylabel(r'$Q_{out}$ ['+reservoir.flux_unit+']')\n", + "ax2.set_xlabel(r'$t$ ['+reservoir.time_unit+']')\n", "ax2.legend()\n", "\n", - "ax3.plot(time_vec[:i_max],pressure_conversion(pressure_vec[:i_max],'Pa',conversion_pressure_unit), label='Pipeline pressure at reservoir')\n", - "ax3.set_ylabel(r'$p_{pipeline}$ ['+conversion_pressure_unit+']')\n", - "ax3.set_xlabel(r'$t$ ['+V.time_unit+']')\n", + "ax3.plot(t_vec,pressure_conversion(pressure_vec,'Pa',pUnit_conv), label='Pipeline pressure at reservoir')\n", + "ax3.set_ylabel(r'$p_{pipeline}$ ['+pUnit_conv+']')\n", + "ax3.set_xlabel(r'$t$ ['+reservoir.time_unit+']')\n", "ax3.legend()\n", "\n", "\n", "fig1.tight_layout() " ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "10.1" - ] - }, - "execution_count": 12, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "V.get_current_level()" - ] } ], "metadata": { diff --git a/Druckrohrleitung/Druckrohrleitung_class_file.py b/Druckrohrleitung/Druckrohrleitung_class_file.py index 1a0ec3d..829ab55 100644 --- a/Druckrohrleitung/Druckrohrleitung_class_file.py +++ b/Druckrohrleitung/Druckrohrleitung_class_file.py @@ -1,5 +1,13 @@ import numpy as np +#importing pressure conversion function +import sys +import os +current = os.path.dirname(os.path.realpath(__file__)) +parent = os.path.dirname(current) +sys.path.append(parent) +from functions.pressure_conversion import pressure_conversion + class Druckrohrleitung_class: # units acceleration_unit = r'$\mathrm{m}/\mathrm{s}^2$' @@ -13,72 +21,66 @@ class Druckrohrleitung_class: velocity_unit = r'$\mathrm{m}/\mathrm{s}$' # for flux and pressure propagation volume_unit = r'$\mathrm{m}^3$' - acceleration_unit_print = 'm/s²' - angle_unit_print = 'rad' - area_unit_print = 'm²' - density_unit_print = 'kg/m³' - flux_unit_print = 'm³/s' - length_unit_print = 'm' - time_unit_print = 's' - velocity_unit_print = 'm/s' # for flux and pressure propagation - volume_unit_print = 'm³' + acceleration_unit_disp = 'm/s²' + angle_unit_disp = 'rad' + area_unit_disp = 'm²' + density_unit_disp = 'kg/m³' + flux_unit_disp = 'm³/s' + length_unit_disp = 'm' + time_unit_disp = 's' + velocity_unit_disp = 'm/s' # for flux and pressure propagation + volume_unit_disp = 'm³' + + g = 9.81 # init - def __init__(self,total_length,diameter,number_segments,pipeline_angle,Darcy_friction_factor,rho=1000,g=9.81): + def __init__(self,total_length,diameter,number_segments,pipeline_angle,Darcy_friction_factor,pw_vel,timestep,pressure_unit_disp,rho=1000): self.length = total_length # total length of the pipeline self.dia = diameter # diameter of the pipeline self.n_seg = number_segments # number of segments for the method of characteristics self.angle = pipeline_angle # angle of the pipeline - self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient + self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient + self.c = pw_vel + self.dt = timestep self.density = rho # density of the liquid in the pipeline - self.g = g # gravitational acceleration self.A = (diameter/2)**2*np.pi self.dx = total_length/number_segments # length of each segment - self.l_vec = np.arange(0,(number_segments+1),1)*self.dx # vector giving the distance from each node to the start of the pipeline + self.x_vec = np.arange(0,(number_segments+1),1)*self.dx # vector giving the distance from each node to the start of the pipeline - # initialize for get_info method - self.c = '--' - self.dt = '--' + self.pressure_unit_disp = pressure_unit_disp # setter - def set_pressure_propagation_velocity(self,c): - self.c = c # propagation velocity of the pressure wave - self.dt = self.dx/c # timestep derived from c, demanded by the method of characteristics - - def set_number_of_timesteps(self,number_timesteps): - self.nt = number_timesteps # number of timesteps - if self.c == '--': - raise Exception('Please set the pressure propagation velocity before setting the number of timesteps.') - else: - self.t_vec = np.arange(0,self.nt*self.dt,self.dt) - def set_initial_pressure(self,pressure): # initialize the pressure distribution in the pipeline if np.size(pressure) == 1: - self.p0 = np.full_like(self.l_vec,pressure) - elif np.size(pressure) == np.size(self.l_vec): - self.p0 = pressure + p0 = np.full_like(self.x_vec,pressure) + elif np.size(pressure) == np.size(self.x_vec): + p0 = pressure else: - raise Exception('Unable to assign initial pressure. Input has to be of size 1 or' + np.size(self.l_vec)) + raise Exception('Unable to assign initial pressure. Input has to be of size 1 or' + np.size(self.x_vec)) #initialize the vectors in which the old and new pressures are stored for the method of characteristics - self.p_old = self.p0.copy() - self.p = self.p0.copy() + self.p_old = p0.copy() + self.p = p0.copy() + self.p_min = p0.copy() + self.p_max = p0.copy() def set_initial_flow_velocity(self,velocity): # initialize the velocity distribution in the pipeline if np.size(velocity) == 1: - self.v0 = np.full_like(self.l_vec,velocity) - elif np.size(velocity) == np.size(self.l_vec): - self.v0 = velocity + v0 = np.full_like(self.x_vec,velocity) + elif np.size(velocity) == np.size(self.x_vec): + v0 = velocity else: - raise Exception('Unable to assign initial velocity. Input has to be of size 1 or' + np.size(self.l_vec)) + raise Exception('Unable to assign initial velocity. Input has to be of size 1 or' + np.size(self.x_vec)) #initialize the vectors in which the old and new velocities are stored for the method of characteristics - self.v_old = self.v0.copy() - self.v = self.v0.copy() + self.v_old = v0.copy() + self.v = v0.copy() + self.v_min = v0.copy() + self.v_max = v0.copy() def set_boundary_conditions_next_timestep(self,p_reservoir,v_turbine): # derived from the method of characteristics, one can set the boundary conditions for the pressures and flow velocities at the reservoir and the turbine @@ -112,16 +114,16 @@ class Druckrohrleitung_class: self.p[0] = p_boundary_res self.p[-1] = p_boundary_tur - def set_steady_state(self,ss_flux,ss_level_reservoir,area_reservoir,pl_vec,h_vec): + def set_steady_state(self,ss_flux,ss_level_reservoir,area_reservoir,x_vec,h_vec): # set the pressure and velocity distributions, that allow a constant flow of water from the (steady-state) reservoir to the (steady-state) turbine # the flow velocity is given by the constant flow through the pipe - ss_v0 = np.full(self.n_seg+1,ss_flux/self.A) + ss_v0 = np.full_like(self.x_vec,ss_flux/self.A) # the static pressure is given by static state pressure of the reservoir, corrected for the hydraulic head of the pipe and friction losses ss_v_in_res = abs(ss_flux/area_reservoir) 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*pl_vec/self.dia*self.density/2*ss_v0**2) + ss_pressure = ss_pressure_res+(self.density*self.g*h_vec)-(self.f_D*x_vec/self.dia*self.density/2*ss_v0**2) self.set_initial_flow_velocity(ss_v0) self.set_initial_pressure(ss_pressure) @@ -135,30 +137,61 @@ class Druckrohrleitung_class: # :<10 pads the self.value to be 10 characters wide print_str = (f"The pipeline has the following attributes: {new_line}" f"----------------------------- {new_line}" - f"Length = {self.length:<10} {self.length_unit_print} {new_line}" - f"Diameter = {self.dia:<10} {self.length_unit_print} {new_line}" + f"Length = {self.length:<10} {self.length_unit_disp} {new_line}" + f"Diameter = {self.dia:<10} {self.length_unit_disp} {new_line}" f"Number of segments = {self.n_seg:<10} {new_line}" f"Number of nodes = {self.n_seg+1:<10} {new_line}" - f"Length per segments = {self.dx:<10} {self.length_unit_print} {new_line}" - f"Pipeline angle = {round(self.angle,3):<10} {self.angle_unit_print} {new_line}" + f"Length per segments = {self.dx:<10} {self.length_unit_disp} {new_line}" + f"Pipeline angle = {round(self.angle,3):<10} {self.angle_unit_disp} {new_line}" f"Pipeline angle = {angle_deg}° {new_line}" f"Darcy friction factor = {self.f_D:<10} {new_line}" - f"Density of liquid = {self.density:<10} {self.density_unit_print} {new_line}" - f"Pressure wave vel. = {self.c:<10} {self.velocity_unit_print} {new_line}" - f"Simulation timestep = {self.dt:<10} {self.time_unit_print} {new_line}" + f"Density of liquid = {self.density:<10} {self.density_unit_disp} {new_line}" + f"Pressure wave vel. = {self.c:<10} {self.velocity_unit_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_print} {new_line}" + f"Total simulation time = {self.nt*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") print(print_str) - def get_current_pressure_distribution(self): - return self.p + def get_current_pressure_distribution(self,disp=False): + if disp == True: + return pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_disp) + elif disp == False: + return self.p def get_current_velocity_distribution(self): return self.v + def get_current_flux_distribution(self): + return self.v*self.A + + def get_lowest_pressure_per_node(self,disp=False): + if disp == True: + return pressure_conversion(self.p_min,self.pressure_unit,self.pressure_unit_disp) + elif disp == False: + return self.p_min + + def get_highest_pressure_per_node(self,disp=False): + if disp == True: + return pressure_conversion(self.p_max,self.pressure_unit,self.pressure_unit_disp) + elif disp == False: + return self.p_max + + def get_lowest_velocity_per_node(self): + return self.v_min + + def get_highest_velocity_per_node(self): + return self.v_max + + def get_lowest_flux_per_node(self): + return self.v_min*self.A + + def get_highest_flux_per_node(self): + return self.v_max*self.A + + def timestep_characteristic_method(self): # 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 @@ -180,6 +213,12 @@ class Druckrohrleitung_class: self.p[i] = 0.5*(self.p_old[i+1]+self.p_old[i-1]) - 0.5*rho*c*(self.v_old[i+1]-self.v_old[i-1]) \ +f_D*rho*c*dt/(4*D)*(abs(self.v_old[i+1])*self.v_old[i+1]-abs(self.v_old[i-1])*self.v_old[i-1]) + # update overall min and max values for pressure and velocity per node + self.p_min = np.minimum(self.p_min,self.p) + self.p_max = np.maximum(self.p_max,self.p) + self.v_min = np.minimum(self.v_min,self.v) + self.v_max = np.maximum(self.v_max,self.v) + # prepare for next call # use .copy() to write data to another memory location and avoid the usual python reference pointer # else one can overwrite data by accidient and change two variables at once without noticing diff --git a/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb b/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb index 83c6b8f..242f710 100644 --- a/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb +++ b/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 9, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -22,102 +22,103 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ - "%matplotlib qt5\n", - "#define constants pipe\n", + "# define constants\n", "\n", - "g = 9.81 # gravitational acceleration [m/s²]\n", - "rho = 1000. # density of water [kg/m³]\n", - "\n", - "L = 1000. # length of pipeline [m]\n", - "D = 0.9 # pipe diameter [m]\n", - "h_res = 10. # water level in upstream reservoir [m]\n", - "n = 50 # number of pipe segments in discretization\n", - "nt = 5000 # number of time steps after initial conditions\n", - "f_D = 0.01 # Darcy friction factor\n", - "c = 400. # propagation velocity of the pressure wave [m/s]\n", - "h_pipe = 105. # hydraulic head without reservoir [m] \n", - "alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n", + " # for physics\n", + "g = 9.81 # [m/s²] gravitational acceleration \n", + "rho = 1000. # [kg/m³] density of water \n", + "pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", + "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n", "\n", "\n", - "# preparing the discretization and initial conditions\n", - "initial_flux = 0.8 # m³/s\n", - "initial_level = h_res # m\n", - "dx = L/n # length of each pipe segment\n", - "dt = dx/c # timestep according to method of characterisitics\n", - "nn = n+1 # number of nodes\n", - "pl_vec = np.arange(0,nn,1)*dx # pl = pipe-length. position of the nodes on the pipeline\n", - "t_vec = np.arange(0,nt,1)*dt # time vector\n", - "h_vec = np.arange(0,nn,1)*h_pipe/n # hydraulic head of pipeline at each node\n", + " # for Turbine\n", + "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", + "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", + "Tur_closingTime = 90. # [s] closing time of turbine\n", "\n", "\n", - "# define constants reservoir\n", - "conversion_pressure_unit = 'mWS'\n", + " # for PI controller\n", + "Con_targetLevel = 8. # [m]\n", + "Con_K_p = 0.1 # [-] proportional constant of PI controller\n", + "Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\n", + "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n", "\n", - "area_base = 75. # m²\n", - "area_pipe = (D/2)**2*np.pi # m²\n", - "critical_level_low = 0. # m\n", - "critical_level_high = 100. # m\n", "\n", - "# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n", - "nt_eRK4 = 1 # number of simulation steps of reservoir in between timesteps of pipeline \n", - "simulation_timestep = dt/nt_eRK4" + " # for pipeline\n", + "Pip_length = (535.+478.) # [m] length of pipeline\n", + "Pip_dia = 0.9 # [m] diameter of pipeline\n", + "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", + "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", + "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", + "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", + "Pip_f_D = 0.014 # [-] Darcy friction factor\n", + "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", + " # derivatives of the pipeline constants\n", + "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", + "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", + "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", + "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", + "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", + "\n", + "\n", + " # for reservoir\n", + "Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n", + "Res_area_out = Pip_area # [m²] outflux area of the reservoir, given by pipeline area\n", + "Res_level_crit_lo = 0. # [m] for yet-to-be-implemented warnings\n", + "Res_level_crit_hi = np.inf # [m] for yet-to-be-implemented warnings\n", + "Res_dt_approx = 1e-3 # [s] approx. timestep of reservoir time evolution to ensure numerical stability (see Res_nt why approx.)\n", + "Res_nt = max(1,int(Pip_dt//Res_dt_approx)) # [1] number of timesteps of the reservoir time evolution within one timestep of the pipeline\n", + "Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n", + "\n", + " # for general simulation\n", + "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", + "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", + "simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", + "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", + "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ - "V = Ausgleichsbecken_class(area_base,area_pipe,critical_level_low,critical_level_high,simulation_timestep)\n", - "V.set_steady_state(initial_flux,initial_level,conversion_pressure_unit)\n", + "# create objects\n", "\n", - "pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n", - "pipe.set_pressure_propagation_velocity(c)\n", - "pipe.set_number_of_timesteps(nt)\n", - "pipe.set_steady_state(initial_flux,initial_level,area_base,pl_vec,h_vec)" + "# Upstream reservoir\n", + "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n", + "reservoir.set_steady_state(flux_init,level_init)\n", + "\n", + "# pipeline\n", + "pipe = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_n_seg,Pip_angle,Pip_f_D,Pip_pw_vel,Pip_dt,pUnit_conv,rho)\n", + "pipe.set_steady_state(flux_init,level_init,Res_area_base,Pip_x_vec,Pip_h_vec)\n" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "The current attributes are: \n", - "----------------------------- \n", - "Current level = 10.0 m\n", - "Volume in reservoir = -- m³ \n", - "Current influx = 0.8 m³/s \n", - "Current outflux = 0.8 m³/s \n", - "Current outflux vel = 1.258 m/s \n", - "Current pipe pressure = 9.844 mWS \n", - "----------------------------- \n", - "\n" - ] - } - ], + "outputs": [], "source": [ - "V.get_info()" + "reservoir.get_info(full=True)\n", + "pipe.get_info(full=True)" ] }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# initialization for timeloop\n", "\n", "level_vec = np.zeros_like(t_vec)\n", - "level_vec[0] = V.get_current_level()\n", + "level_vec[0] = reservoir.get_current_level()\n", "\n", "# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n", "v_old = pipe.get_current_velocity_distribution()\n", @@ -141,21 +142,23 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ + "%matplotlib qt5\n", "fig1,axs1 = plt.subplots(2,1)\n", "axs1[0].set_title('Pressure distribution in pipeline')\n", "axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs1[0].set_ylabel(r'$p$ [mWS]')\n", - "lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,'Pa',conversion_pressure_unit),marker='.')\n", - "axs1[0].set_ylim([0.9*np.min(pressure_conversion(p_old,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_old,'Pa',conversion_pressure_unit))])\n", + "axs1[0].set_ylim([0.9*np.min(pressure_conversion(p_old,'Pa',pUnit_conv)),1.1*np.max(pressure_conversion(p_old,'Pa',pUnit_conv))])\n", + "lo_00, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,'Pa',pUnit_conv),marker='.')\n", "\n", "axs1[1].set_title('Velocity distribution in pipeline')\n", "axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs1[1].set_ylabel(r'$v$ [m/s]')\n", - "lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n", + "lo_01, = axs1[1].plot(Pip_x_vec,v_old,marker='.')\n", + "axs1[1].autoscale()\n", "# axs1[1].set_ylim([0.9*np.min(v_old),1.1*np.max(v_boundary_res)])\n", "\n", "fig1.tight_layout()\n", @@ -164,25 +167,24 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ - "\n", - "for it_pipe in range(1,nt):\n", + "for it_pipe in range(1,nt+1):\n", "# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n", " # set initial conditions for the reservoir time evolution calculted with e-RK4\n", - " V.set_pressure(p_old[0])\n", - " V.set_outflux(v_old[0]*area_pipe)\n", + " reservoir.set_pressure(p_old[0],display_warning=False)\n", + " reservoir.set_outflux(v_old[0]*Pip_area,display_warning=False)\n", " # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n", - " for it_res in range(nt_eRK4):\n", - " V.timestep_reservoir_evolution() \n", - " level_vec[it_pipe] = V.get_current_level() \n", + " for it_res in range(Res_nt):\n", + " reservoir.timestep_reservoir_evolution() \n", + " level_vec[it_pipe] = reservoir.get_current_level() \n", "\n", " \n", " # set boundary conditions for the next timestep of the characteristic method\n", - " p_boundary_res[it_pipe] = V.get_current_pressure()\n", - " v_boundary_tur[it_pipe] = initial_flux/area_pipe\n", + " p_boundary_res[it_pipe] = reservoir.get_current_pressure()\n", + " v_boundary_tur[it_pipe] = flux_init/Pip_area\n", "\n", " # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n", " pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n", @@ -202,28 +204,30 @@ " lo_01.remove()\n", " # lo_02.remove()\n", " # plot new pressure and velocity distribution in the pipeline\n", - " lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,'Pa', conversion_pressure_unit),marker='.',c='blue')\n", - " lo_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n", + " lo_00, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,'Pa', pUnit_conv),marker='.',c='blue')\n", + " lo_01, = axs1[1].plot(Pip_x_vec,v_old,marker='.',c='blue')\n", " \n", " fig1.suptitle(str(round(t_vec[it_pipe],2)) + '/' + str(round(t_vec[-1],2)))\n", " fig1.canvas.draw()\n", " fig1.tight_layout()\n", " plt.pause(0.000001)\n", - "\n" + "\n", + "reservoir.get_info(full=True)\n", + "pipe.get_info(full=True)" ] }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "fig2,axs2 = plt.subplots(2,2)\n", "axs2[0,0].set_title('Pressure Reservoir')\n", - "axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit))\n", + "axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,pUnit_calc,pUnit_conv))\n", "axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2[0,0].set_ylabel(r'$p$ [mWS]')\n", - "axs2[0,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit))])\n", + "axs2[0,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_res,pUnit_calc,pUnit_conv)),1.1*np.max(pressure_conversion(p_boundary_res,pUnit_calc,pUnit_conv))])\n", "\n", "axs2[0,1].set_title('Velocity Reservoir')\n", "axs2[0,1].plot(t_vec,v_boundary_res)\n", @@ -232,16 +236,16 @@ "axs2[0,1].set_ylim([0.9*np.min(v_boundary_res),1.1*np.max(v_boundary_res)])\n", "\n", "axs2[1,0].set_title('Pressure Turbine')\n", - "axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit))\n", + "axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,pUnit_calc,pUnit_conv))\n", "axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2[1,0].set_ylabel(r'$p$ [mWS]')\n", - "axs2[1,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit))])\n", + "axs2[1,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_tur,pUnit_calc,pUnit_conv)),1.1*np.max(pressure_conversion(p_boundary_tur,pUnit_calc,pUnit_conv))])\n", "\n", "axs2[1,1].set_title('Velocity Turbine')\n", "axs2[1,1].plot(t_vec,v_boundary_tur)\n", "axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", - "axs2[1,1].set_ylim([0.9*np.min(v_boundary_tur),1.1*np.max(v_boundary_tur)])\n", + "axs2[1,1].set_ylim([0.95*np.min(v_boundary_tur),1.05*np.max(v_boundary_tur)])\n", "\n", "fig2.tight_layout()\n", "plt.show()" diff --git a/Regler/Pegelregler_test.ipynb b/Regler/Pegelregler_test.ipynb index 081cd94..4d16b95 100644 --- a/Regler/Pegelregler_test.ipynb +++ b/Regler/Pegelregler_test.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 27, "metadata": {}, "outputs": [], "source": [ @@ -23,85 +23,108 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 28, "metadata": {}, "outputs": [], "source": [ - "#define constants\n", + "# define constants\n", "\n", - "#Turbine\n", - "Q_nenn = 0.85 # m³/s\n", - "p_nenn = pressure_conversion(10.6,'bar','Pa')\n", - "closing_time = 480. #s\n", + " # for physics\n", + "g = 9.81 # [m/s²] gravitational acceleration \n", + "rho = 1000. # [kg/m³] density of water \n", + "pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", + "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n", "\n", - "# physics\n", - "g = 9.81 # gravitational acceleration [m/s²]\n", - "rho = 1000. # density of water [kg/m³]\n", "\n", - "# define controller constants\n", - "target_level = 8. # m\n", - "Kp = 0.01\n", - "Ti = 3600.\n", - "deadband_range = 0.05 # m\n", + " # for Turbine\n", + "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", + "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", + "Tur_closingTime = 90. # [s] closing time of turbine\n", "\n", - "# reservoir\n", - "initial_level = target_level\n", - "initial_influx = Q_nenn/2 # initial influx of volume to the reservoir [m³/s]\n", - "initial_pressure_unit = 'Pa' # DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", - "conversion_pressure_unit = 'bar' # for pressure conversion in print statements and plot labels\n", - "area_base = 74. # total base are of the cuboid reservoir [m²] \n", - "area_outflux = 1. # outflux area of the reservoir, given by pipeline area [m²]\n", - "critical_level_low = 0. # for yet-to-be-implemented warnings[m]\n", - "critical_level_high = np.inf # for yet-to-be-implemented warnings[m]\n", "\n", - "p0 = rho*g*initial_level-0.5*rho*(initial_influx/area_outflux)**2\n", + " # for PI controller\n", + "Con_targetLevel = 8. # [m]\n", + "Con_K_p = 0.1 # [-] proportional constant of PI controller\n", + "Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\n", + "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n", "\n", - "# offset the pressure in front of the turbine to get realisitc fluxes\n", - "h_fict = 100\n", - "offset_pressure = rho*g*h_fict\n", "\n", - "t_max = 1e4 #s\n", - "dt = 1e-2 # simulation timestep\n", - "nt = int(t_max//dt) # number of simulation steps of reservoir in between timesteps of pipeline \n", + " # for pipeline\n", + "Pip_length = (535.+478.) # [m] length of pipeline\n", + "Pip_dia = 0.9 # [m] diameter of pipeline\n", + "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", + "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", + "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", + "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", + "Pip_f_D = 0.014 # [-] Darcy friction factor\n", + "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", + " # derivatives of the pipeline constants\n", + "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", + "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", + "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", + "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", + "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", "\n", - "t_vec = np.arange(0,nt+1,1)*dt\n", - "\n" + "\n", + " # for reservoir\n", + "Res_area_base = 10. # [m²] total base are of the cuboid reservoir \n", + "Res_area_out = Pip_area # [m²] outflux area of the reservoir, given by pipeline area\n", + "Res_level_crit_lo = 0. # [m] for yet-to-be-implemented warnings\n", + "Res_level_crit_hi = np.inf # [m] for yet-to-be-implemented warnings\n", + "Res_dt_approx = 1e-3 # [s] approx. timestep of reservoir time evolution to ensure numerical stability (see Res_nt why approx.)\n", + "Res_nt = max(1,int(Pip_dt//Res_dt_approx)) # [1] number of timesteps of the reservoir time evolution within one timestep of the pipeline\n", + "Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n", + "\n", + " # for general simulation\n", + "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", + "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", + "simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", + "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", + "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "# create objects\n", + "offset_pressure = pressure_conversion(Pip_head,'mws',pUnit_calc)\n", "\n", - "V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,dt)\n", - "V.set_steady_state(initial_influx,initial_level,conversion_pressure_unit)\n", + "# Upstream reservoir\n", + "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n", + "reservoir.set_steady_state(flux_init,level_init)\n", "\n", - "T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time,dt)\n", - "T1.set_steady_state(initial_influx,p0+offset_pressure)\n", + "# downstream turbine\n", + "turbine = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime,Pip_dt,pUnit_conv)\n", + "turbine.set_steady_state(flux_init,reservoir.get_current_pressure()+offset_pressure)\n", "\n", - "Pegelregler = PI_controller_class(target_level,deadband_range,Kp,Ti,dt)" + "\n", + "# level controll\n", + "level_control = PI_controller_class(Con_targetLevel,Con_deadbandRange,Con_K_p,Con_T_i,Pip_dt)\n", + "level_control.set_control_variable(turbine.get_current_LA(),display_warning=False)\n" ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 30, "metadata": {}, "outputs": [], "source": [ - "level_vec = np.full(nt+1,V.level)\n", - "LA_ist_vec = np.full(nt+1,T1.LA)\n", - "LA_soll_vec = np.full(nt+1,T1.LA)\n", - "Q_vec = np.full(nt+1,initial_influx)\n", - "\n", - "Pegelregler.control_variable = T1.get_current_LA()" + "level_vec = np.zeros_like(t_vec)\n", + "level_vec[0] = level_init\n", + "LA_ist_vec = np.zeros_like(t_vec)\n", + "LA_ist_vec[0] = turbine.get_current_LA()\n", + "LA_soll_vec = np.zeros_like(t_vec)\n", + "LA_soll_vec[0] = turbine.get_current_LA()\n", + "Q_vec = np.zeros_like(t_vec)\n", + "Q_vec[0] = turbine.get_current_Q()" ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 31, "metadata": {}, "outputs": [ { @@ -109,105 +132,20 @@ "output_type": "stream", "text": [ "0.0\n", - "100.0\n", - "200.0\n", - "300.0\n", - "400.0\n", - "500.0\n", - "600.0\n", - "700.0\n", - "800.0\n", - "900.0\n", - "1000.0\n", - "1100.0\n", - "1200.0\n", - "1300.0\n", - "1400.0\n", - "1500.0\n", - "1600.0\n", - "1700.0\n", - "1800.0\n", - "1900.0\n", - "2000.0\n", - "2100.0\n", - "2200.0\n", - "2300.0\n", - "2400.0\n", - "2500.0\n", - "2600.0\n", - "2700.0\n", - "2800.0\n", - "2900.0\n", - "3000.0\n", - "3100.0\n", - "3200.0\n", - "3300.0\n", - "3400.0\n", - "3500.0\n", - "3600.0\n", - "3700.0\n", - "3800.0\n", - "3900.0\n", - "4000.0\n", - "4100.0\n", - "4200.0\n", - "4300.0\n", - "4400.0\n", - "4500.0\n", - "4600.0\n", - "4700.0\n", - "4800.0\n", - "4900.0\n", - "5000.0\n", - "5100.0\n", - "5200.0\n", - "5300.0\n", - "5400.0\n", - "5500.0\n", - "5600.0\n", - "5700.0\n", - "5800.0\n", - "5900.0\n", - "6000.0\n", - "6100.0\n", - "6200.0\n", - "6300.0\n", - "6400.0\n", - "6500.0\n", - "6600.0\n", - "6700.0\n", - "6800.0\n", - "6900.0\n", - "7000.0\n", - "7100.0\n", - "7200.0\n", - "7300.0\n", - "7400.0\n", - "7500.0\n", - "7600.0\n", - "7700.0\n", - "7800.0\n", - "7900.0\n", - "8000.0\n", - "8100.0\n", - "8200.0\n", - "8300.0\n", - "8400.0\n", - "8500.0\n", - "8600.0\n", - "8700.0\n", - "8800.0\n", - "8900.0\n", - "9000.0\n", - "9100.0\n", - "9200.0\n", - "9300.0\n", - "9400.0\n", - "9500.0\n", - "9600.0\n", - "9700.0\n", - "9800.0\n", - "9900.0\n" + "40.52\n", + "81.04\n", + "121.56\n", + "162.08\n", + "202.6\n", + "243.12\n", + "283.64\n", + "324.16\n", + "364.68\n", + "405.2\n", + "445.72\n", + "486.24\n", + "526.76\n", + "567.28\n" ] } ], @@ -216,34 +154,34 @@ "\n", "for i in range(nt+1):\n", "\n", - " if np.mod(i,1e4) == 0:\n", + " if np.mod(i,1e3) == 0:\n", " print(t_vec[i])\n", "\n", - " if i == 0.4*(nt+1):\n", - " V.set_influx(0.)\n", + " if i > 0.1*(nt+1):\n", + " reservoir.set_influx(0.)\n", "\n", - " p = V.get_current_pressure()\n", - " Pegelregler.update_control_variable(V.level)\n", - " LA_soll = Pegelregler.get_current_control_variable()\n", - " T1.update_LA(LA_soll)\n", - " T1.set_pressure(p+offset_pressure)\n", + " p = reservoir.get_current_pressure()\n", + " level_control.update_control_variable(reservoir.level)\n", + " LA_soll = level_control.get_current_control_variable()\n", + " turbine.update_LA(LA_soll)\n", + " turbine.set_pressure(p+offset_pressure)\n", " LA_soll_vec[i] = LA_soll\n", - " LA_ist_vec[i] = T1.get_current_LA()\n", - " Q_vec[i] = T1.get_current_Q()\n", + " LA_ist_vec[i] = turbine.get_current_LA()\n", + " Q_vec[i] = turbine.get_current_Q()\n", "\n", " \n", - " V.set_outflux(Q_vec[i])\n", + " reservoir.set_outflux(Q_vec[i],display_warning=False)\n", "\n", - " V.timestep_reservoir_evolution() \n", - " \n", - " level_vec[i] = V.get_current_level()\n", + " for it_res in range(Res_nt):\n", + " reservoir.timestep_reservoir_evolution() \n", + " level_vec[i] = reservoir.get_current_level()\n", " \n", " " ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 32, "metadata": {}, "outputs": [], "source": [ @@ -256,12 +194,12 @@ "axs1[0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs1[0].set_ylabel(r'$h$ [$\\mathrm{m}$]')\n", "axs1[0].plot(t_vec,level_vec)\n", - "axs1[0].set_ylim([0*initial_level,1.5*initial_level])\n", + "axs1[0].set_ylim([0*level_init,1.5*level_init])\n", "axs1[1].set_title('Flux')\n", "axs1[1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs1[1].set_ylabel(r'$Q$ [$\\mathrm{m} / \\mathrm{s}^3$]')\n", "axs1[1].plot(t_vec,Q_vec)\n", - "axs1[1].set_ylim([0,2*initial_influx])\n", + "axs1[1].set_ylim([0,2*flux_init])\n", "axs1[2].set_title('LA')\n", "axs1[2].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs1[2].set_ylabel(r'$LA$ [%]')\n", @@ -271,27 +209,6 @@ "fig1.tight_layout()\n", "fig1.show()\n" ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[]" - ] - }, - "execution_count": 7, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "fig2 = plt.figure()\n", - "plt.plot(t_vec,Pegelregler.get_error_history())" - ] } ], "metadata": { diff --git a/Regler/Regler_class_file.py b/Regler/Regler_class_file.py index dcbf58a..88c298b 100644 --- a/Regler/Regler_class_file.py +++ b/Regler/Regler_class_file.py @@ -84,17 +84,17 @@ class PI_controller_class: # use a list to be able to append more easily - will get converted to np.array when needed self.error_history = [0] - self.control_variable = -99 - self.cv_lower_limit = lower_limit # limits for the controll variable self.cv_upper_limit = upper_limit # limits for the controll variable # setter + + def set_setpoint(self,setpoint): self.SP = setpoint def set_control_variable(self,control_variable, display_warning=True): - if display_warning == True and self.control_variable != -99: + if display_warning == True: print('WARNING! You are setting the control variable of the PI controller manually \ and are not using the .update_controll_variable() method') self.control_variable = control_variable diff --git a/Turbinen/Turbinen_class_file.py b/Turbinen/Turbinen_class_file.py index 30e37a0..21f895e 100644 --- a/Turbinen/Turbinen_class_file.py +++ b/Turbinen/Turbinen_class_file.py @@ -1,8 +1,10 @@ -from time import time import numpy as np + #importing pressure conversion function import sys import os + +from pyparsing import alphanums current = os.path.dirname(os.path.realpath(__file__)) parent = os.path.dirname(current) sys.path.append(parent) @@ -11,7 +13,7 @@ from functions.pressure_conversion import pressure_conversion class Francis_Turbine: # units # make sure that units and print units are the same - # units are used to label graphs and print units are used to have a bearable format when using pythons print() + # 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$' flux_unit = r'$\mathrm{m}^3/\mathrm{s}$' LA_unit = '%' @@ -20,30 +22,28 @@ class Francis_Turbine: velocity_unit = r'$\mathrm{m}/\mathrm{s}$' volume_unit = r'$\mathrm{m}^3$' - density_unit_print = 'kg/m³' - flux_unit_print = 'm³/s' - LA_unit_print = '%' - pressure_unit_print = 'mWS' - time_unit_print = 's' - velocity_unit_print = 'm/s' - volume_unit_print = 'm³' + density_unit_disp = 'kg/m³' + flux_unit_disp = 'm³/s' + LA_unit_disp = '%' + time_unit_disp = 's' + velocity_unit_disp = 'm/s' + volume_unit_disp = 'm³' g = 9.81 # m/s² gravitational acceleration # init - def __init__(self, Q_nenn,p_nenn,t_closing=-1.,timestep=-1.): + def __init__(self, Q_nenn,p_nenn,t_closing,timestep,pressure_unit_disp): self.Q_n = Q_nenn # nominal flux self.p_n = p_nenn # nominal pressure self.LA_n = 1. # 100% # nominal Leitapparatöffnung - h = pressure_conversion(p_nenn,'Pa','MWs') # nominal pressure in terms of hydraulic head - self.A = Q_nenn/(np.sqrt(2*self.g*h)*0.98) # Ersatzfläche - self.dt = timestep # simulation timestep - self.t_c = t_closing # closing time + self.t_c = t_closing # closing time self.d_LA_max_dt = 1/t_closing # maximal change of LA per second + self.pressure_unit_disp = pressure_unit_disp + # initialize for get_info() - parameters will be converted to display -1 if not overwritten - self.p = pressure_conversion(-1,self.pressure_unit_print,self.pressure_unit) + self.p = pressure_conversion(-1,self.pressure_unit_disp,self.pressure_unit) self.Q = -1. self.LA = -0.01 @@ -54,19 +54,22 @@ class Francis_Turbine: self.LA = LA # warn user, that the .set_LA() method should not be used ot set LA manually if display_warning == True: - print('Consider using the .update_LA() method instead of setting LA manually') - - def set_timestep(self,timestep,display_warning=True): - # set Leitapparatöffnung - self.dt = time - # warn user, that the .set_LA() method should not be used ot set LA manually - if display_warning == True: - print('WARNING: You are changing the timestep of the turbine simulation. This has implications on the simulated closing speed!') + print('You are setting the guide vane opening of the turbine manually. \n \ + This is not an intended use of this method. \n \ + Refer to the .update_LA() method instead.') def set_pressure(self,pressure): # set pressure in front of the turbine self.p = pressure + def set_steady_state(self,ss_flux,ss_pressure): + # calculate and set steady state LA, that allows the flow of ss_flux at ss_pressure through the + # turbine at the steady state LA + ss_LA = self.LA_n*ss_flux/self.Q_n*np.sqrt(self.p_n/ss_pressure) + if ss_LA < 0 or ss_LA > 1: + raise Exception('LA out of range [0;1]') + self.set_LA(ss_LA,display_warning=False) + #getter def get_current_Q(self): # return the flux through the turbine, based on the current pressure in front @@ -80,10 +83,13 @@ class Francis_Turbine: def get_current_LA(self): return self.LA + def get_current_pressure(self): + return pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_disp) + def get_info(self, full = False): new_line = '\n' - p = pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_print) - p_n = pressure_conversion(self.p_n,self.pressure_unit,self.pressure_unit_print) + p = pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_disp) + p_n = pressure_conversion(self.p_n,self.pressure_unit,self.pressure_unit_disp) if full == True: @@ -91,33 +97,34 @@ class Francis_Turbine: print_str = (f"Turbine has the following attributes: {new_line}" f"----------------------------- {new_line}" f"Type = Francis {new_line}" - f"Nominal flux = {self.Q_n:<10} {self.flux_unit_print} {new_line}" - f"Nominal pressure = {round(p_n,3):<10} {self.pressure_unit_print}{new_line}" - f"Nominal LA = {self.LA_n*100:<10} {self.LA_unit_print} {new_line}" - f"Closing time = {self.t_c:<10} {self.time_unit_print} {new_line}" - f"Current flux = {self.Q:<10} {self.flux_unit_print} {new_line}" - f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}" - f"Current LA = {self.LA*100:<10} {self.LA_unit_print} {new_line}" - f"Simulation timestep = {self.dt:<10} {self.time_unit_print} {new_line}" + f"Nominal flux = {self.Q_n:<10} {self.flux_unit_disp} {new_line}" + f"Nominal pressure = {round(p_n,3):<10} {self.pressure_unit_disp}{new_line}" + f"Nominal LA = {self.LA_n*100:<10} {self.LA_unit_disp} {new_line}" + f"Closing time = {self.t_c:<10} {self.time_unit_disp} {new_line}" + f"Current flux = {self.Q:<10} {self.flux_unit_disp} {new_line}" + f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_disp} {new_line}" + f"Current LA = {self.LA*100:<10} {self.LA_unit_disp} {new_line}" + f"Simulation timestep = {self.dt:<10} {self.time_unit_disp} {new_line}" f"----------------------------- {new_line}") else: # :<10 pads the self.value to be 10 characters wide print_str = (f"The current attributes are: {new_line}" f"----------------------------- {new_line}" - f"Current flux = {self.Q:<10} {self.flux_unit_print} {new_line}" - f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}" - f"Current LA = {self.LA*100:<10} {self.LA_unit_print} {new_line}" + f"Current flux = {self.Q:<10} {self.flux_unit_disp} {new_line}" + f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_disp} {new_line}" + f"Current LA = {self.LA*100:<10} {self.LA_unit_disp} {new_line}" f"----------------------------- {new_line}") print(print_str) -# methods +# update methods def update_LA(self,LA_soll): # update the Leitappartöffnung and consider the restrictions of the closing time of the turbine - LA_diff = self.LA-LA_soll # calculate the difference to the target LA - LA_diff_max = self.d_LA_max_dt*self.dt # calculate the maximum change in LA based on the given timestep + LA_diff = self.LA-LA_soll # calculate the difference to the target LA + LA_diff_max = self.d_LA_max_dt*self.dt # calculate the maximum possible change in LA based on the given timestep LA_diff = np.sign(LA_diff)*np.min(np.abs([LA_diff,LA_diff_max])) # calulate the correct change in LA + # make sure that the LA is not out of the range [0;1] LA_new = self.LA-LA_diff if LA_new < 0.: LA_new = 0. @@ -125,10 +132,42 @@ class Francis_Turbine: LA_new = 1. self.set_LA(LA_new,display_warning=False) - def set_steady_state(self,ss_flux,ss_pressure): - # calculate and set steady state LA, that allows the flow of ss_flux at ss_pressure through the - # turbine at the steady state LA - ss_LA = self.LA_n*ss_flux/self.Q_n*np.sqrt(self.p_n/ss_pressure) - if ss_LA < 0 or ss_LA > 1: - raise Exception('LA out of range [0;1]') - self.set_LA(ss_LA,display_warning=False) +# methods + 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 compatible. + eps = 1e-12 # convergence criterion: iteration change < eps + iteration_change = 1. # change in Q from one iteration to the next + i = 0 # safety variable. break loop if it exceeds 1e6 iterations + g = self.g # gravitational acceleration + p = convergence_parameters[0] # pressure at second to last node (see method of characterisctics - boundary condidtions) + v = convergence_parameters[1] # velocity at second to last node (see method of characterisctics - boundary condidtions) + D = convergence_parameters[2] # diameter of the pipeline + area_pipe = convergence_parameters[3] # area of the pipeline + alpha = convergence_parameters[4] # elevation angle of the pipeline + f_D = convergence_parameters[5] # Darcy friction coefficient + c = convergence_parameters[6] # pressure wave propagtation velocity + rho = convergence_parameters[7] # density of the liquid + dt = convergence_parameters[8] # timestep of the characteristic method + + p_old = self.get_current_pressure() + Q_old = self.get_current_Q() + v_old = Q_old/area_pipe + + + while iteration_change > eps: + self.set_pressure(p_old) + Q_new = self.get_current_Q() + v_new = Q_new/area_pipe + p_new = p-rho*c*(v_old-v)+rho*c*dt*g*np.sin(alpha)-f_D*rho*c*dt/(2*D)*abs(v)*v + + iteration_change = abs(Q_old-Q_new) + Q_old = Q_new.copy() + p_old = p_new.copy() + v_old = v_new.copy() + i = i+1 + if i == 1e6: + print('did not converge') + break + # print(i) + self.Q = Q_new \ No newline at end of file diff --git a/Turbinen/Turbinen_test_steady_state.ipynb b/Turbinen/Turbinen_test_steady_state.ipynb new file mode 100644 index 0000000..4037691 --- /dev/null +++ b/Turbinen/Turbinen_test_steady_state.ipynb @@ -0,0 +1,370 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from Turbinen_class_file import Francis_Turbine\n", + "\n", + "import sys\n", + "import os\n", + "current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n", + "parent = os.path.dirname(current)\n", + "sys.path.append(parent)\n", + "from functions.pressure_conversion import pressure_conversion\n", + "from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class\n", + "from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class\n", + "from Regler.Regler_class_file import PI_controller_class" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "# define constants\n", + "\n", + " # for physics\n", + "g = 9.81 # [m/s²] gravitational acceleration \n", + "rho = 1000. # [kg/m³] density of water \n", + "pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", + "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n", + "\n", + "\n", + " # for Turbine\n", + "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", + "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", + "Tur_closingTime = 90. # [s] closing time of turbine\n", + "\n", + "\n", + " # for PI controller\n", + "Con_targetLevel = 8. # [m]\n", + "Con_K_p = 0.1 # [-] proportional constant of PI controller\n", + "Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\n", + "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n", + "\n", + "\n", + " # for pipeline\n", + "Pip_length = (535.+478.) # [m] length of pipeline\n", + "Pip_dia = 0.9 # [m] diameter of pipeline\n", + "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", + "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", + "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", + "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", + "Pip_f_D = 0.014 # [-] Darcy friction factor\n", + "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", + " # derivatives of the pipeline constants\n", + "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", + "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", + "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", + "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", + "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", + "\n", + "\n", + " # for reservoir\n", + "Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n", + "Res_area_out = Pip_area # [m²] outflux area of the reservoir, given by pipeline area\n", + "Res_level_crit_lo = 0. # [m] for yet-to-be-implemented warnings\n", + "Res_level_crit_hi = np.inf # [m] for yet-to-be-implemented warnings\n", + "Res_dt_approx = 1e-3 # [s] approx. timestep of reservoir time evolution to ensure numerical stability (see Res_nt why approx.)\n", + "Res_nt = max(1,int(Pip_dt//Res_dt_approx)) # [1] number of timesteps of the reservoir time evolution within one timestep of the pipeline\n", + "Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n", + "\n", + " # for general simulation\n", + "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", + "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", + "simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", + "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", + "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "# create objects\n", + "\n", + "# Upstream reservoir\n", + "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n", + "reservoir.set_steady_state(flux_init,level_init)\n", + "\n", + "# pipeline\n", + "pipe = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_n_seg,Pip_angle,Pip_f_D,Pip_pw_vel,Pip_dt,pUnit_conv,rho)\n", + "pipe.set_steady_state(flux_init,level_init,Res_area_base,Pip_x_vec,Pip_h_vec)\n", + "\n", + "# downstream turbine\n", + "turbine = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime,Pip_dt,pUnit_conv)\n", + "turbine.set_steady_state(flux_init,pipe.get_current_pressure_distribution()[-1])\n", + "\n", + "# influx setting turbine\n", + "turbine_in = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime/2,Pip_dt,pUnit_conv)\n", + "turbine_in.set_steady_state(flux_init,Tur_p_nenn)\n", + "\n", + "# level controll\n", + "level_control = PI_controller_class(Con_targetLevel,Con_deadbandRange,Con_K_p,Con_T_i,Pip_dt)\n", + "level_control.set_control_variable(turbine.get_current_LA(),display_warning=False)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "# initialization for Timeloop\n", + "\n", + "v_old = pipe.get_current_velocity_distribution()\n", + "v_min = pipe.get_current_velocity_distribution()\n", + "v_max = pipe.get_current_velocity_distribution()\n", + "Q_old = pipe.get_current_flux_distribution()\n", + "Q_min = pipe.get_current_flux_distribution()\n", + "Q_max = pipe.get_current_flux_distribution()\n", + "p_old = pipe.get_current_pressure_distribution()\n", + "p_min = pipe.get_current_pressure_distribution()\n", + "p_max = pipe.get_current_pressure_distribution()\n", + "\n", + "Q_in_vec = np.zeros_like(t_vec)\n", + "Q_in_vec[0] = flux_init\n", + "\n", + "v_boundary_res = np.zeros_like(t_vec)\n", + "v_boundary_tur = np.zeros_like(t_vec)\n", + "Q_boundary_res = np.zeros_like(t_vec)\n", + "Q_boundary_tur = np.zeros_like(t_vec)\n", + "p_boundary_res = np.zeros_like(t_vec)\n", + "p_boundary_tur = np.zeros_like(t_vec)\n", + "\n", + "level_vec = np.full_like(t_vec,level_init) # level at the end of each pipeline timestep\n", + "volume_vec = np.full_like(t_vec,reservoir.get_current_volume()) # volume at the end of each pipeline timestep\n", + "\n", + "v_boundary_res[0] = v_old[0]\n", + "v_boundary_tur[0] = v_old[-1] \n", + "Q_boundary_res[0] = Q_old[0]\n", + "Q_boundary_tur[0] = Q_old[-1]\n", + "p_boundary_res[0] = p_old[0]\n", + "p_boundary_tur[0] = p_old[-1]\n", + "\n", + "LA_soll_vec = np.full_like(t_vec,turbine.get_current_LA())\n", + "LA_ist_vec = np.full_like(t_vec,turbine.get_current_LA())\n", + "\n", + "LA_soll_vec2 = np.full_like(t_vec,turbine_in.get_current_LA())\n" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib qt5\n", + "# Con_T_ime loop\n", + "\n", + "# create a figure and subplots to display the velocity and pressure distribution across the pipeline in each pipeline step\n", + "fig1,axs1 = plt.subplots(2,1)\n", + "fig1.suptitle(str(0) +' s / '+str(round(t_vec[-1],2)) + ' s' )\n", + "axs1[0].set_title('Pressure distribution in pipeline')\n", + "axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", + "axs1[0].set_ylabel(r'$p$ ['+pUnit_conv+']')\n", + "axs1[1].set_title('Flux distribution in pipeline')\n", + "axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", + "axs1[1].set_ylabel(r'$Q$ [$\\mathrm{m}^3 / \\mathrm{s}$]')\n", + "lo_p, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,pUnit_calc, pUnit_conv),marker='.')\n", + "lo_q, = axs1[1].plot(Pip_x_vec,Q_old,marker='.')\n", + "lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n", + "lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", + "lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n", + "lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n", + "\n", + "axs1[0].autoscale()\n", + "axs1[1].autoscale()\n", + "\n", + "fig1.tight_layout()\n", + "fig1.show()\n", + "plt.pause(1)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "convergence_parameters = [p_old[-2],v_old[-2],Pip_dia,Pip_area,Pip_angle,Pip_f_D,Pip_pw_vel,rho,Pip_dt]\n", + "\n", + "# loop through Con_T_ime steps of the pipeline\n", + "for it_pipe in range(1,nt+1):\n", + "\n", + " turbine_in.update_LA(LA_soll_vec2[it_pipe])\n", + " turbine_in.set_pressure(Tur_p_nenn)\n", + " Q_in_vec[it_pipe] = turbine_in.get_current_Q()\n", + " reservoir.set_influx(Q_in_vec[it_pipe])\n", + "\n", + "# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n", + " # set initial condition for the reservoir Con_T_ime evolution calculted with e-RK4\n", + " reservoir.set_pressure(p_old[0],display_warning=False)\n", + " reservoir.set_outflux(Q_old[0],display_warning=False)\n", + " # calculate the Con_T_ime evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n", + " for it_res in range(Res_nt):\n", + " reservoir.timestep_reservoir_evolution() \n", + " level_vec[it_pipe] = reservoir.get_current_level() \n", + " volume_vec[it_pipe] = reservoir.get_current_volume() \n", + "\n", + " # get the control variable\n", + " level_control.update_control_variable(level_vec[it_pipe])\n", + " LA_soll_vec[it_pipe] = level_control.get_current_control_variable()\n", + " \n", + " # change the Leitapparatöffnung based on the target value\n", + " turbine.update_LA(LA_soll_vec[it_pipe])\n", + " LA_ist_vec[it_pipe] = turbine.get_current_LA()\n", + "\n", + " # set boundary condition for the next timestep of the characterisCon_T_ic method\n", + " turbine.set_pressure(p_old[-1])\n", + " convergence_parameters[0] = p_old[-2]\n", + " convergence_parameters[1] = v_old[-2]\n", + " turbine.converge(convergence_parameters)\n", + " p_boundary_res[it_pipe] = reservoir.get_current_pressure()\n", + " v_boundary_tur[it_pipe] = 1/Pip_area*turbine.get_current_Q()\n", + " Q_boundary_tur[it_pipe] = turbine.get_current_Q()\n", + "\n", + " # the the boundary condition in the pipe.object and thereby calculate boundary pressure at turbine\n", + " pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n", + " pipe.v[0] = (0.8*pipe.v[0]+0.2*reservoir.get_current_outflux()/Res_area_out)\n", + " p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\n", + " v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n", + " Q_boundary_res[it_pipe] = pipe.get_current_flux_distribution()[0]\n", + "\n", + " # perform the next timestep via the characterisCon_T_ic method\n", + " pipe.timestep_characteristic_method()\n", + "\n", + " # prepare for next loop\n", + " p_old = pipe.get_current_pressure_distribution()\n", + " v_old = pipe.get_current_velocity_distribution()\n", + " Q_old = pipe.get_current_flux_distribution()\n", + "\n", + "\n", + " # plot some stuff\n", + " # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n", + " lo_p.remove()\n", + " lo_pmin.remove()\n", + " lo_pmax.remove()\n", + " lo_q.remove()\n", + " lo_qmin.remove()\n", + " lo_qmax.remove()\n", + " # plot new pressure and velocity distribution in the pipeline\n", + " lo_p, = axs1[0].plot(Pip_x_vec,pipe.get_current_pressure_distribution(disp=True),marker='.',c='blue')\n", + " lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n", + " lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", + " lo_q, = axs1[1].plot(Pip_x_vec,pipe.get_current_flux_distribution(),marker='.',c='blue')\n", + " lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n", + " lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n", + " fig1.suptitle(str(round(t_vec[it_pipe],2))+ ' s / '+str(round(t_vec[-1],2)) + ' s' )\n", + " fig1.canvas.draw()\n", + " fig1.tight_layout()\n", + " fig1.show()\n", + " plt.pause(0.001) " + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "# plot Con_T_ime evolution of boundary pressure and velocity as well as the reservoir level\n", + "\n", + "fig2,axs2 = plt.subplots(1,1)\n", + "axs2.set_title('Level and Volume reservoir')\n", + "axs2.plot(t_vec,level_vec,label='level')\n", + "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", + "axs2.set_ylabel(r'$h$ [m]')\n", + "x_twin_00 = axs2.twinx()\n", + "x_twin_00.set_ylabel(r'$V$ [$\\mathrm{m}^3$]')\n", + "x_twin_00.plot(t_vec,volume_vec)\n", + "axs2.legend()\n", + "\n", + "fig2,axs2 = plt.subplots(1,1)\n", + "axs2.set_title('LA')\n", + "axs2.plot(t_vec,100*LA_soll_vec,label='Target')\n", + "axs2.plot(t_vec,100*LA_ist_vec,label='Actual')\n", + "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", + "axs2.set_ylabel(r'$LA$ [%]')\n", + "axs2.legend()\n", + "\n", + "fig2,axs2 = plt.subplots(1,1)\n", + "axs2.set_title('Pressure reservoir and turbine')\n", + "axs2.plot(t_vec,pressure_conversion(p_boundary_res,pUnit_calc, pUnit_conv),label='Reservoir')\n", + "axs2.plot(t_vec,pressure_conversion(p_boundary_tur,pUnit_calc, pUnit_conv),label='Turbine')\n", + "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", + "axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n", + "axs2.legend()\n", + "\n", + "fig2,axs2 = plt.subplots(1,1)\n", + "axs2.set_title('Fluxes')\n", + "axs2.plot(t_vec,Q_boundary_res,label='Outflux')\n", + "axs2.plot(t_vec,Q_in_vec,label='Influx')\n", + "axs2.plot(t_vec,Q_boundary_tur,label='Flux Turbine')\n", + "axs2.set_ylim(-2*Tur_Q_nenn,+2*Tur_Q_nenn)\n", + "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", + "axs2.set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n", + "axs2.legend()\n", + "\n", + "fig2,axs2 = plt.subplots(1,1)\n", + "axs2.set_title('Min and Max Pressure')\n", + "axs2.plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n", + "axs2.plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", + "axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", + "axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n", + "\n", + "fig2,axs2 = plt.subplots(1,1)\n", + "axs2.set_title('Min and Max Fluxes')\n", + "axs2.plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n", + "axs2.plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n", + "axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", + "axs2.set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n", + "\n", + "# axs2[0,1].legend()\n", + "# axs2[1,0].legend()\n", + "# axs2[1,1].legend()\n", + "# # axs2[2,0].legend()\n", + "# # axs2[2,1].legend()\n", + "\n", + "\n", + "fig2.tight_layout()\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3.8.13 ('Georg_DT_Slot3')", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.13" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "84fb123bdc47ab647d3782661abcbe80fbb79236dd2f8adf4cef30e8755eb2cd" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/Turbinen/convergence_turbine.py b/Turbinen/old/convergence_turbine.py similarity index 99% rename from Turbinen/convergence_turbine.py rename to Turbinen/old/convergence_turbine.py index 5930c3e..6945a6d 100644 --- a/Turbinen/convergence_turbine.py +++ b/Turbinen/old/convergence_turbine.py @@ -1,4 +1,3 @@ -from time import time import numpy as np #importing pressure conversion function import sys @@ -59,7 +58,7 @@ class Francis_Turbine_test: def set_timestep(self,timestep,display_warning=True): # set Leitapparatöffnung - self.dt = time + self.dt = timestep # warn user, that the .set_LA() method should not be used ot set LA manually if display_warning == True: print('WARNING: You are changing the timestep of the turbine simulation. This has implications on the simulated closing speed!') diff --git a/Turbinen/turbine_convergence_test.ipynb b/Turbinen/old/turbine_convergence_test.ipynb similarity index 75% rename from Turbinen/turbine_convergence_test.ipynb rename to Turbinen/old/turbine_convergence_test.ipynb index 5e8c347..2f71293 100644 --- a/Turbinen/turbine_convergence_test.ipynb +++ b/Turbinen/old/turbine_convergence_test.ipynb @@ -18,7 +18,8 @@ "sys.path.append(parent)\n", "from functions.pressure_conversion import pressure_conversion\n", "from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class\n", - "from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class" + "from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class\n", + "from Regler.Regler_class_file import PI_controller_class" ] }, { @@ -32,26 +33,32 @@ "#Turbine\n", "Q_nenn = 0.85 # m³/s\n", "p_nenn = pressure_conversion(10.6,'bar','Pa')\n", - "closing_time = 5 #s\n", + "closing_time = 90. #s\n", "\n", "#define constants pipe\n", "\n", "g = 9.81 # gravitational acceleration [m/s²]\n", "rho = 1000. # density of water [kg/m³]\n", "\n", - "L = 1000. # length of pipeline [m]\n", + "# define controller constants\n", + "target_level = 8. # m\n", + "Kp = 0.1\n", + "Ti = 1000.\n", + "deadband_range = 0.05 # m\n", + "\n", + "L = 535.+478. # length of pipeline [m]\n", "D = 0.9 # pipe diameter [m]\n", - "h_res = 10. # water level in upstream reservoir [m]\n", + "h_res = target_level # water level in upstream reservoir [m]\n", "n = 50 # number of pipe segments in discretization\n", "nt = 10000 # number of time steps after initial conditions\n", - "f_D = 0.01 # Darcy friction factor\n", + "f_D = 0.014 # Darcy friction factor\n", "c = 400. # propagation velocity of the pressure wave [m/s]\n", "h_pipe = 105. # hydraulic head without reservoir [m] \n", "alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n", "\n", "\n", "# preparing the discretization and initial conditions\n", - "initial_flux = 0.8 # m³/s\n", + "initial_flux = Q_nenn/1.1 # m³/s\n", "initial_level = h_res # m\n", "dx = L/n # length of each pipe segment\n", "dt = dx/c # timestep according to method of characterisitics\n", @@ -70,7 +77,7 @@ "critical_level_high = 100. # m\n", "\n", "# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n", - "nt_eRK4 = 1 # number of simulation steps of reservoir in between timesteps of pipeline \n", + "nt_eRK4 = 100 # number of simulation steps of reservoir in between timesteps of pipeline \n", "simulation_timestep = dt/nt_eRK4" ] }, @@ -91,7 +98,13 @@ "\n", "initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n", "T1 = Francis_Turbine_test(Q_nenn,p_nenn,closing_time,timestep=dt)\n", - "T1.set_steady_state(initial_flux,initial_pressure_turbine)" + "T1.set_steady_state(initial_flux,initial_pressure_turbine)\n", + "\n", + "T_in = Francis_Turbine_test(Q_nenn,p_nenn,closing_time/2,timestep=dt)\n", + "T_in.set_steady_state(initial_flux,p_nenn)\n", + "\n", + "Pegelregler = PI_controller_class(target_level,deadband_range,Kp,Ti,dt)\n", + "Pegelregler.control_variable = T1.get_current_LA()" ] }, { @@ -122,7 +135,16 @@ "v_boundary_res[0] = v_old[0]\n", "v_boundary_tur[0] = v_old[-1] \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", + "LA_soll_vec = np.full_like(t_vec,T1.get_current_LA())\n", + "LA_ist_vec = np.full_like(t_vec,T1.get_current_LA())\n", + "\n", + "LA_soll_vec2 = np.full_like(t_vec,T_in.get_current_LA())\n", + "LA_soll_vec2[500:1000] = 0.\n", + "LA_soll_vec2[1000:1500] = 1. \n", + "LA_soll_vec2[1500:2000] = 0.\n", + "LA_soll_vec2[2000:2500] = 0.5 " ] }, { @@ -136,13 +158,14 @@ "axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs1[0].set_ylabel(r'$p$ [mWS]')\n", "lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,'Pa',conversion_pressure_unit),marker='.')\n", - "axs1[0].set_ylim([0.9*np.min(pressure_conversion(p_old,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_old,'Pa',conversion_pressure_unit))])\n", "\n", "axs1[1].set_title('Velocity distribution in pipeline')\n", "axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs1[1].set_ylabel(r'$v$ [m/s]')\n", "lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n", - "# axs1[1].set_ylim([0.9*np.min(v_old),1.1*np.max(v_boundary_res)])\n", + "\n", + "axs1[0].autoscale()\n", + "axs1[1].autoscale()\n", "\n", "fig1.tight_layout()\n", "plt.pause(1)" @@ -157,6 +180,11 @@ "\n", "for it_pipe in range(1,nt):\n", "# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n", + " \n", + " T_in.update_LA(LA_soll_vec2[it_pipe])\n", + " T_in.set_pressure(p_nenn)\n", + " V.set_influx(T_in.get_current_Q())\n", + "\n", " # set initial conditions for the reservoir time evolution calculted with e-RK4\n", " V.set_pressure(p_old[0])\n", " V.set_outflux(v_old[0]*area_pipe)\n", @@ -165,16 +193,23 @@ " V.timestep_reservoir_evolution() \n", " level_vec[it_pipe] = V.get_current_level() \n", "\n", + " # get the control variable\n", + " Pegelregler.update_control_variable(level_vec[it_pipe])\n", + " LA_soll_vec[it_pipe] = Pegelregler.get_current_control_variable()\n", " \n", + " # change the Leitapparatöffnung based on the target value\n", + " T1.update_LA(LA_soll_vec[it_pipe])\n", + " LA_ist_vec[it_pipe] = T1.get_current_LA()\n", + "\n", " # set boundary conditions for the next timestep of the characteristic method\n", - " p_boundary_res[it_pipe] = V.get_current_pressure()\n", " T1.set_pressure(p_old[-1])\n", " T1.converge(area_pipe,p_old[-2],v_old[-2],alpha,f_D,dt)\n", + " p_boundary_res[it_pipe] = V.get_current_pressure()\n", " v_boundary_tur[it_pipe] = T1.get_current_Q()/area_pipe\n", "\n", " # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n", " pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n", - " pipe.v[0] = (pipe.v[0]+V.get_current_outflux()/area_pipe)/2\n", + " pipe.v[0] = (0.8*pipe.v[0]+0.2*V.get_current_outflux()/area_pipe)\n", " p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\n", " v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n", "\n", @@ -202,58 +237,47 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ - "fig2,axs2 = plt.subplots(2,2)\n", - "axs2[0,0].set_title('Pressure Reservoir')\n", - "axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit))\n", - "axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[0,0].set_ylabel(r'$p$ [mWS]')\n", - "axs2[0,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit))])\n", + "# plot time evolution of boundary pressure and velocity as well as the reservoir level\n", "\n", - "axs2[0,1].set_title('Velocity Reservoir')\n", + "fig2,axs2 = plt.subplots(3,2)\n", + "axs2[0,0].set_title('Pressure reservoir')\n", + "axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,'Pa', conversion_pressure_unit))\n", + "axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", + "axs2[0,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", + "\n", + "axs2[0,1].set_title('Velocity reservoir')\n", "axs2[0,1].plot(t_vec,v_boundary_res)\n", + "axs2[0,1].set_ylim(-2*Q_nenn,+2*Q_nenn)\n", "axs2[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", - "axs2[0,1].set_ylim([0.9*np.min(v_boundary_res),1.1*np.max(v_boundary_res)])\n", "\n", - "axs2[1,0].set_title('Pressure Turbine')\n", - "axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit))\n", + "axs2[1,0].set_title('Pressure turbine')\n", + "axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,'Pa', conversion_pressure_unit))\n", "axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[1,0].set_ylabel(r'$p$ [mWS]')\n", - "axs2[1,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit))])\n", + "axs2[1,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", "\n", - "axs2[1,1].set_title('Velocity Turbine')\n", + "axs2[1,1].set_title('Velocity turbine')\n", "axs2[1,1].plot(t_vec,v_boundary_tur)\n", "axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", - "axs2[1,1].set_ylim([0.9*np.min(v_boundary_tur),1.1*np.max(v_boundary_tur)])\n", "\n", + "axs2[2,0].set_title('Level reservoir')\n", + "axs2[2,0].plot(t_vec,level_vec)\n", + "axs2[2,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", + "axs2[2,0].set_ylabel(r'$h$ [m]')\n", + "\n", + "axs2[2,1].set_title('LA')\n", + "axs2[2,1].plot(t_vec,100*LA_soll_vec)\n", + "axs2[2,1].plot(t_vec,100*LA_ist_vec)\n", + "axs2[2,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", + "axs2[2,1].set_ylabel(r'$LA$ [%]')\n", "fig2.tight_layout()\n", "plt.show()" ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[]" - ] - }, - "execution_count": 8, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "plt.plot(level_vec)" - ] } ], "metadata": { diff --git a/Untertweng.ipynb b/Untertweng.ipynb deleted file mode 100644 index bbd030a..0000000 --- a/Untertweng.ipynb +++ /dev/null @@ -1,360 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "\n", - "from functions.pressure_conversion import pressure_conversion\n", - "from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class\n", - "from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class\n", - "from Turbinen.Turbinen_class_file import Francis_Turbine" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "#define constants\n", - "\n", - "#Turbine\n", - "Q_nenn = 0.85 # m³/s\n", - "p_nenn = pressure_conversion(10.6,'bar','Pa')\n", - "closing_time = 5 #s\n", - "\n", - "# physics\n", - "g = 9.81 # gravitational acceleration [m/s²]\n", - "rho = 1000. # density of water [kg/m³]\n", - "\n", - "# pipeline\n", - "L = 535.+478. # length of pipeline [m]\n", - "D = 0.9 # pipe diameter [m]\n", - "A_pipe = D**2/4*np.pi # pipeline area\n", - "h_pipe = 105 # hydraulic head without reservoir [m] \n", - "alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n", - "n = 50 # number of pipe segments in discretization # initial flow velocity [m/s]\n", - "f_D = 0.014 # Darcy friction factor\n", - "c = 500. # propagation velocity of the pressure wave [m/s]\n", - "# consider prescribing a total simulation time and deducting the number of timesteps from that\n", - "nt = 2500 # number of time steps after initial conditions\n", - "\n", - "# derivatives of the pipeline constants\n", - "dx = L/n # length of each pipe segment\n", - "dt = dx/c # timestep according to method of characterisitics\n", - "nn = n+1 # number of nodes\n", - "initial_level = 8. # water level in upstream reservoir [m]\n", - "pl_vec = np.arange(0,nn,1)*dx # pl = pipe-length. position of the nodes on the pipeline\n", - "t_vec = np.arange(0,nt+1)*dt # time vector\n", - "h_vec = np.arange(0,nn,1)*h_pipe/n # hydraulic head of pipeline at each node \n", - "\n", - "\n", - "\n", - "# reservoir\n", - "# replace influx by vector\n", - "initial_flux = Q_nenn/1.1 # initial influx of volume to the reservoir [m³/s]\n", - "initial_pressure_unit = 'Pa' # DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", - "conversion_pressure_unit = 'bar' # for pressure conversion in print statements and plot labels\n", - "area_base = 74. # total base are of the cuboid reservoir [m²] \n", - "area_outflux = A_pipe # outlfux area of the reservoir, given by pipeline area [m²]\n", - "critical_level_low = 0. # for yet-to-be-implemented warnings[m]\n", - "critical_level_high = np.inf # for yet-to-be-implemented warnings[m]\n", - "\n", - "# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n", - "nt_eRK4 = 100 # number of simulation steps of reservoir in between timesteps of pipeline \n", - "simulation_timestep = dt/nt_eRK4\n", - "\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "# create objects\n", - "\n", - "V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,simulation_timestep)\n", - "V.set_steady_state(initial_flux,initial_level,conversion_pressure_unit)\n", - "\n", - "\n", - "pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n", - "pipe.set_pressure_propagation_velocity(c)\n", - "pipe.set_number_of_timesteps(nt)\n", - "pipe.set_steady_state(initial_flux,initial_level,area_base,pl_vec,h_vec)\n", - "\n", - "initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n", - "\n", - "T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time,timestep=dt)\n", - "T1.set_steady_state(initial_flux,initial_pressure_turbine)\n" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "# initialization for timeloop\n", - "\n", - "# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n", - "v_old = pipe.get_current_velocity_distribution()\n", - "p_old = pipe.get_current_pressure_distribution()\n", - "\n", - "# prepare the vectors in which the temporal evolution of the boundary conditions are stored\n", - " # keep in mind, that the velocity at the turbine and the pressure at the reservoir follow from boundary conditions\n", - " # reservoir level and flow through turbine\n", - " # the pressure at the turbine and the velocity at the reservoir are calculated from the method of characteristics\n", - "v_boundary_res = np.zeros_like(t_vec)\n", - "v_boundary_tur = np.zeros_like(t_vec)\n", - "p_boundary_res = np.zeros_like(t_vec)\n", - "p_boundary_tur = np.zeros_like(t_vec)\n", - "\n", - "# prepare the vectors that store the temporal evolution of the level in the reservoir\n", - "level_vec = np.full(nt+1,initial_level) # level at the end of each pipeline timestep\n", - "\n", - "# set the boundary conditions for the first timestep\n", - "v_boundary_res[0] = v_old[0]\n", - "v_boundary_tur[0] = v_old[-1] \n", - "p_boundary_res[0] = p_old[0]\n", - "p_boundary_tur[0] = p_old[-1]\n", - "\n", - "LA_soll_vec = np.full_like(t_vec,T1.LA)\n", - "LA_soll_vec[500:]= 0\n", - "LA_ist_vec = np.full_like(t_vec,T1.LA)\n", - "\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib qt5\n", - "# time loop\n", - "\n", - "# create a figure and subplots to display the velocity and pressure distribution across the pipeline in each pipeline step\n", - "fig1,axs1 = plt.subplots(2,1)\n", - "fig1.suptitle(str(0) +' s / '+str(round(t_vec[-1],2)) + ' s' )\n", - "axs1[0].set_title('Pressure distribution in pipeline')\n", - "axs1[1].set_title('Velocity distribution in pipeline')\n", - "axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", - "axs1[0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", - "axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", - "axs1[1].set_ylabel(r'$v$ [$\\mathrm{m} / \\mathrm{s}$]')\n", - "lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.')\n", - "lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n", - "axs1[0].autoscale()\n", - "axs1[1].autoscale()\n", - "\n", - "fig1.tight_layout()\n", - "fig1.show()\n", - "plt.pause(1)\n" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "# loop through time steps of the pipeline\n", - "for it_pipe in range(1,pipe.nt+1):\n", - "\n", - " if it_pipe == 250:\n", - " V.set_influx(0.)\n", - "\n", - "# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n", - " # set initial conditions for the reservoir time evolution calculted with e-RK4\n", - " V.set_pressure(p_old[0])\n", - " V.set_outflux(v_old[0]*area_outflux)\n", - " # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n", - " for it_res in range(nt_eRK4):\n", - " V.timestep_reservoir_evolution() \n", - " level_vec[it_pipe] = V.get_current_level() \n", - "\n", - " # change the Leitapparatöffnung based on the target value\n", - " T1.update_LA(LA_soll_vec[it_pipe])\n", - " T1.set_pressure(p_old[-1])\n", - "\n", - " LA_ist_vec[it_pipe] = T1.LA\n", - "\n", - " # set boundary conditions for the next timestep of the characteristic method\n", - " p_boundary_res[it_pipe] = V.get_current_pressure()\n", - " v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_current_Q()\n", - "\n", - " # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n", - " pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n", - " p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\n", - " v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n", - "\n", - " # perform the next timestep via the characteristic method\n", - " pipe.timestep_characteristic_method()\n", - "\n", - " # prepare for next loop\n", - " p_old = pipe.get_current_pressure_distribution()\n", - " v_old = pipe.get_current_velocity_distribution()\n", - "\n", - " # plot some stuff\n", - " # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n", - " lo_00.remove()\n", - " lo_01.remove()\n", - " # lo_02.remove()\n", - " # plot new pressure and velocity distribution in the pipeline\n", - " lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.',c='blue')\n", - " lo_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n", - " # lo_02, = axs1[2].plot(level_vec_2,c='blue')\n", - " fig1.suptitle(str(round(t_vec[it_pipe],2))+ ' s / '+str(round(t_vec[-1],2)) + ' s' )\n", - " fig1.canvas.draw()\n", - " fig1.tight_layout()\n", - " fig1.show()\n", - " plt.pause(0.001) \n", - "\n", - " \n", - " " - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "0.0\n" - ] - } - ], - "source": [ - "print(V.get_current_influx())" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [], - "source": [ - "# plot time evolution of boundary pressure and velocity as well as the reservoir level\n", - "\n", - "fig2,axs2 = plt.subplots(3,2)\n", - "axs2[0,0].set_title('Pressure reservoir')\n", - "axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,initial_pressure_unit, conversion_pressure_unit))\n", - "axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[0,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", - "\n", - "axs2[0,1].set_title('Velocity reservoir')\n", - "axs2[0,1].plot(t_vec,v_boundary_res)\n", - "axs2[0,1].set_ylim(-2*Q_nenn,+2*Q_nenn)\n", - "axs2[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", - "\n", - "axs2[1,0].set_title('Pressure turbine')\n", - "axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,initial_pressure_unit, conversion_pressure_unit))\n", - "axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[1,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", - "\n", - "axs2[1,1].set_title('Velocity turbine')\n", - "axs2[1,1].plot(t_vec,v_boundary_tur)\n", - "axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", - "\n", - "axs2[2,0].set_title('Level reservoir')\n", - "axs2[2,0].plot(t_vec,level_vec)\n", - "axs2[2,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[2,0].set_ylabel(r'$h$ [m]')\n", - "\n", - "axs2[2,1].set_title('LA')\n", - "axs2[2,1].plot(t_vec,100*LA_soll_vec)\n", - "axs2[2,1].plot(t_vec,100*LA_ist_vec)\n", - "axs2[2,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[2,1].set_ylabel(r'$LA$ [%]')\n", - "fig2.tight_layout()\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "The cuboid reservoir has the following attributes: \n", - "----------------------------- \n", - "Base area = 74.0 m² \n", - "Outflux area = 0.636 m² \n", - "Current level = 7.875725956447418 m\n", - "Critical level low = 0.0 m \n", - "Critical level high = inf m \n", - "Volume in reservoir = -- m³ \n", - "Current influx = 0.0 m³/s \n", - "Current outflux = -0.1415386124341686 m³/s \n", - "Current outflux vel = -0.222 m/s \n", - "Current pipe pressure = 0.772 bar \n", - "Simulation timestep = 0.0004052 s \n", - "Density of liquid = 1000 kg/m³ \n", - "----------------------------- \n", - "\n", - "9.22707730779877\n", - "10.57842865915012\n", - "11.92978001050147\n", - "13.281131361852822\n", - "14.632482713204173\n", - "15.983834064555523\n", - "17.335185415906874\n", - "18.686536767258225\n", - "20.037888118609576\n", - "21.389239469960927\n" - ] - } - ], - "source": [ - "V.get_info(full=True)\n", - "V.set_outflux(-10.)\n", - "for i in range(10):\n", - " V.level = V.update_level(10.)\n", - " print(V.get_current_level())" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3.8.13 ('DT_Slot_3')", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.13" - }, - "orig_nbformat": 4, - "vscode": { - "interpreter": { - "hash": "4a28055eb8a3160fa4c7e4fca69770c4e0a1add985300856aa3fcf4ce32a2c48" - } - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/Untertweng_mit_Pegelregler.ipynb b/Untertweng_mit_Pegelregler.ipynb index 835e3da..8ea99dd 100644 --- a/Untertweng_mit_Pegelregler.ipynb +++ b/Untertweng_mit_Pegelregler.ipynb @@ -18,67 +18,64 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ - "#define constants\n", + "# define constants\n", "\n", - "#Turbine\n", - "Q_nenn = 0.85 # m³/s\n", - "p_nenn = pressure_conversion(10.6,'bar','Pa')\n", - "closing_time = 30. #s\n", - "\n", - "# physics\n", - "g = 9.81 # gravitational acceleration [m/s²]\n", - "rho = 1000. # density of water [kg/m³]\n", + " # for physics\n", + "g = 9.81 # [m/s²] gravitational acceleration \n", + "rho = 1000. # [kg/m³] density of water \n", + "pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", + "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n", "\n", "\n", - "# define controller constants\n", - "target_level = 8. # m\n", - "Kp = 0.1\n", - "Ti = 7.\n", - "deadband_range = 0.05 # m\n", + " # for Turbine\n", + "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", + "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", + "Tur_closingTime = 90. # [s] closing time of turbine\n", "\n", "\n", - "# pipeline\n", - "L = (535.+478.) # length of pipeline [m]\n", - "D = 0.9 # pipe diameter [m]\n", - "A_pipe = D**2/4*np.pi # pipeline area\n", - "h_pipe = 105 # hydraulic head without reservoir [m] \n", - "alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n", - "n = 50 # number of pipe segments in discretization\n", - "f_D = 0.014 # Darcy friction factor\n", - "c = 500. # propagation velocity of the pressure wave [m/s]\n", - "# consider prescribing a total simulation time and deducting the number of timesteps from that\n", - "nt = 9000 # number of time steps after initial conditions\n", - "\n", - "# derivatives of the pipeline constants\n", - "dx = L/n # length of each pipe segment\n", - "dt = dx/c # timestep according to method of characterisitics\n", - "nn = n+1 # number of nodes\n", - "initial_level = target_level # water level in upstream reservoir [m]\n", - "pl_vec = np.arange(0,nn,1)*dx # pl = pipe-length. position of the nodes on the pipeline\n", - "t_vec = np.arange(0,nt+1)*dt # time vector\n", - "h_vec = np.arange(0,nn,1)*h_pipe/n # hydraulic head of pipeline at each node \n", + " # for PI controller\n", + "Con_targetLevel = 8. # [m]\n", + "Con_K_p = 0.1 # [-] proportional constant of PI controller\n", + "Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\n", + "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n", "\n", "\n", + " # for pipeline\n", + "Pip_length = (535.+478.) # [m] length of pipeline\n", + "Pip_dia = 0.9 # [m] diameter of pipeline\n", + "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", + "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", + "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", + "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", + "Pip_f_D = 0.014 # [-] Darcy friction factor\n", + "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", + " # derivatives of the pipeline constants\n", + "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", + "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", + "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", + "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", + "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", "\n", - "# reservoir\n", - "# replace influx by vector\n", - "initial_flux = Q_nenn/1.1 # initial influx of volume to the reservoir [m³/s]\n", - "initial_pressure_unit = 'Pa' # DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", - "conversion_pressure_unit = 'bar' # for pressure conversion in print statements and plot labels\n", - "area_base = 74. # total base are of the cuboid reservoir [m²] \n", - "area_outflux = A_pipe # outlfux area of the reservoir, given by pipeline area [m²]\n", - "critical_level_low = 0. # for yet-to-be-implemented warnings[m]\n", - "critical_level_high = np.inf # for yet-to-be-implemented warnings[m]\n", "\n", - "# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n", - "nt_eRK4 = 100 # number of simulation steps of reservoir in between timesteps of pipeline \n", - "simulation_timestep = dt/nt_eRK4\n", + " # for reservoir\n", + "Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n", + "Res_area_out = Pip_area # [m²] outflux area of the reservoir, given by pipeline area\n", + "Res_level_crit_lo = 0. # [m] for yet-to-be-implemented warnings\n", + "Res_level_crit_hi = np.inf # [m] for yet-to-be-implemented warnings\n", + "Res_dt_approx = 1e-3 # [s] approx. timestep of reservoir time evolution to ensure numerical stability (see Res_nt why approx.)\n", + "Res_nt = max(1,int(Pip_dt//Res_dt_approx)) # [1] number of timesteps of the reservoir time evolution within one timestep of the pipeline\n", + "Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n", "\n", - "\n" + " # for general simulation\n", + "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", + "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", + "simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", + "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", + "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" ] }, { @@ -89,65 +86,25 @@ "source": [ "# create objects\n", "\n", - "V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,simulation_timestep)\n", - "V.set_steady_state(initial_flux,initial_level,conversion_pressure_unit)\n", + "# Upstream reservoir\n", + "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n", + "reservoir.set_steady_state(flux_init,level_init)\n", "\n", + "# pipeline\n", + "pipe = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_n_seg,Pip_angle,Pip_f_D,Pip_pw_vel,Pip_dt,pUnit_conv,rho)\n", + "pipe.set_steady_state(flux_init,level_init,Res_area_base,Pip_x_vec,Pip_h_vec)\n", "\n", - "pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n", - "pipe.set_pressure_propagation_velocity(c)\n", - "pipe.set_number_of_timesteps(nt)\n", - "pipe.set_steady_state(initial_flux,initial_level,area_base,pl_vec,h_vec)\n", + "# downstream turbine\n", + "turbine = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime,Pip_dt,pUnit_conv)\n", + "turbine.set_steady_state(flux_init,pipe.get_current_pressure_distribution()[-1])\n", "\n", - "initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n", + "# influx setting turbine\n", + "turbine_in = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime/2,Pip_dt,pUnit_conv)\n", + "turbine_in.set_steady_state(flux_init,Tur_p_nenn)\n", "\n", - "T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time,timestep=dt)\n", - "T1.set_steady_state(initial_flux,initial_pressure_turbine)\n", - "\n", - "T_in = Francis_Turbine(Q_nenn,p_nenn,closing_time/2,timestep=dt)\n", - "T_in.set_steady_state(initial_flux,p_nenn)\n", - "\n", - "Pegelregler = PI_controller_class(target_level,deadband_range,Kp,Ti,dt)\n", - "Pegelregler.control_variable = T1.get_current_LA()\n" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "# initialization for timeloop\n", - "\n", - "# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n", - "v_old = pipe.get_current_velocity_distribution()\n", - "p_old = pipe.get_current_pressure_distribution()\n", - "\n", - "# prepare the vectors in which the temporal evolution of the boundary conditions are stored\n", - " # keep in mind, that the velocity at the turbine and the pressure at the reservoir follow from boundary conditions\n", - " # reservoir level and flow through turbine\n", - " # the pressure at the turbine and the velocity at the reservoir are calculated from the method of characteristics\n", - "v_boundary_res = np.zeros_like(t_vec)\n", - "v_boundary_tur = np.zeros_like(t_vec)\n", - "p_boundary_res = np.zeros_like(t_vec)\n", - "p_boundary_tur = np.zeros_like(t_vec)\n", - "\n", - "# prepare the vectors that store the temporal evolution of the level in the reservoir\n", - "level_vec = np.full(nt+1,initial_level) # level at the end of each pipeline timestep\n", - "\n", - "# set the boundary conditions for the first timestep\n", - "v_boundary_res[0] = v_old[0]\n", - "v_boundary_tur[0] = v_old[-1] \n", - "p_boundary_res[0] = p_old[0]\n", - "p_boundary_tur[0] = p_old[-1]\n", - "\n", - "LA_soll_vec = np.full_like(t_vec,T1.get_current_LA())\n", - "LA_ist_vec = np.full_like(t_vec,T1.get_current_LA())\n", - "\n", - "LA_soll_vec2 = np.full_like(t_vec,T_in.get_current_LA())\n", - "LA_soll_vec2[500:1000] = 0.\n", - "LA_soll_vec2[1000:1500] = 1. \n", - "LA_soll_vec2[1500:2000] = 0.\n", - "LA_soll_vec2[2000:2500] = 0.5 \n" + "# level controll\n", + "level_control = PI_controller_class(Con_targetLevel,Con_deadbandRange,Con_K_p,Con_T_i,Pip_dt)\n", + "level_control.set_control_variable(turbine.get_current_LA(),display_warning=False)\n" ] }, { @@ -155,21 +112,74 @@ "execution_count": 5, "metadata": {}, "outputs": [], + "source": [ + "# initialization for Timeloop\n", + "\n", + "v_old = pipe.get_current_velocity_distribution()\n", + "v_min = pipe.get_current_velocity_distribution()\n", + "v_max = pipe.get_current_velocity_distribution()\n", + "Q_old = pipe.get_current_flux_distribution()\n", + "Q_min = pipe.get_current_flux_distribution()\n", + "Q_max = pipe.get_current_flux_distribution()\n", + "p_old = pipe.get_current_pressure_distribution()\n", + "p_min = pipe.get_current_pressure_distribution()\n", + "p_max = pipe.get_current_pressure_distribution()\n", + "\n", + "Q_in_vec = np.zeros_like(t_vec)\n", + "Q_in_vec[0] = flux_init\n", + "\n", + "v_boundary_res = np.zeros_like(t_vec)\n", + "v_boundary_tur = np.zeros_like(t_vec)\n", + "Q_boundary_res = np.zeros_like(t_vec)\n", + "Q_boundary_tur = np.zeros_like(t_vec)\n", + "p_boundary_res = np.zeros_like(t_vec)\n", + "p_boundary_tur = np.zeros_like(t_vec)\n", + "\n", + "level_vec = np.full_like(t_vec,level_init) # level at the end of each pipeline timestep\n", + "volume_vec = np.full_like(t_vec,reservoir.get_current_volume()) # volume at the end of each pipeline timestep\n", + "\n", + "v_boundary_res[0] = v_old[0]\n", + "v_boundary_tur[0] = v_old[-1] \n", + "Q_boundary_res[0] = Q_old[0]\n", + "Q_boundary_tur[0] = Q_old[-1]\n", + "p_boundary_res[0] = p_old[0]\n", + "p_boundary_tur[0] = p_old[-1]\n", + "\n", + "LA_soll_vec = np.full_like(t_vec,turbine.get_current_LA())\n", + "LA_ist_vec = np.full_like(t_vec,turbine.get_current_LA())\n", + "\n", + "LA_soll_vec2 = np.full_like(t_vec,turbine_in.get_current_LA())\n", + "LA_soll_vec2[500:1000] = 0.\n", + "LA_soll_vec2[1000:1500] = 1. \n", + "LA_soll_vec2[1500:2000] = 0.\n", + "LA_soll_vec2[2000:2500] = 0.5 \n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], "source": [ "%matplotlib qt5\n", - "# time loop\n", + "# Con_T_ime loop\n", "\n", "# create a figure and subplots to display the velocity and pressure distribution across the pipeline in each pipeline step\n", "fig1,axs1 = plt.subplots(2,1)\n", "fig1.suptitle(str(0) +' s / '+str(round(t_vec[-1],2)) + ' s' )\n", "axs1[0].set_title('Pressure distribution in pipeline')\n", - "axs1[1].set_title('Velocity distribution in pipeline')\n", "axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", - "axs1[0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", + "axs1[0].set_ylabel(r'$p$ ['+pUnit_conv+']')\n", + "axs1[1].set_title('Flux distribution in pipeline')\n", "axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", - "axs1[1].set_ylabel(r'$v$ [$\\mathrm{m} / \\mathrm{s}$]')\n", - "lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.')\n", - "lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n", + "axs1[1].set_ylabel(r'$Q$ [$\\mathrm{m}^3 / \\mathrm{s}$]')\n", + "lo_p, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,pUnit_calc, pUnit_conv),marker='.')\n", + "lo_q, = axs1[1].plot(Pip_x_vec,Q_old,marker='.')\n", + "lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n", + "lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", + "lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n", + "lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n", + "\n", "axs1[0].autoscale()\n", "axs1[1].autoscale()\n", "\n", @@ -180,136 +190,153 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ - "error_vec = np.zeros_like(t_vec)\n", - "# loop through time steps of the pipeline\n", - "for it_pipe in range(1,pipe.nt+1):\n", + "convergence_parameters = [p_old[-2],v_old[-2],Pip_dia,Pip_area,Pip_angle,Pip_f_D,Pip_pw_vel,rho,Pip_dt]\n", "\n", - " T_in.update_LA(LA_soll_vec2[it_pipe])\n", - " T_in.set_pressure(p_nenn)\n", - " V.set_influx(T_in.get_current_Q())\n", + "# loop through Con_T_ime steps of the pipeline\n", + "for it_pipe in range(1,nt+1):\n", + "\n", + " turbine_in.update_LA(LA_soll_vec2[it_pipe])\n", + " turbine_in.set_pressure(Tur_p_nenn)\n", + " Q_in_vec[it_pipe] = turbine_in.get_current_Q()\n", + " reservoir.set_influx(Q_in_vec[it_pipe])\n", "\n", "# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n", - " # set initial conditions for the reservoir time evolution calculted with e-RK4\n", - " V.set_pressure(p_old[0])\n", - " V.set_outflux(v_old[0]*area_outflux)\n", - " # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n", - " for it_res in range(nt_eRK4):\n", - " V.timestep_reservoir_evolution() \n", - " level_vec[it_pipe] = V.get_current_level() \n", + " # set initial condition for the reservoir Con_T_ime evolution calculted with e-RK4\n", + " reservoir.set_pressure(p_old[0],display_warning=False)\n", + " reservoir.set_outflux(Q_old[0],display_warning=False)\n", + " # calculate the Con_T_ime evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n", + " for it_res in range(Res_nt):\n", + " reservoir.timestep_reservoir_evolution() \n", + " level_vec[it_pipe] = reservoir.get_current_level() \n", + " volume_vec[it_pipe] = reservoir.get_current_volume() \n", "\n", " # get the control variable\n", - " Pegelregler.update_control_variable(level_vec[it_pipe])\n", - " LA_soll_vec[it_pipe] = Pegelregler.get_current_control_variable()\n", + " level_control.update_control_variable(level_vec[it_pipe])\n", + " LA_soll_vec[it_pipe] = level_control.get_current_control_variable()\n", " \n", " # change the Leitapparatöffnung based on the target value\n", - " T1.update_LA(LA_soll_vec[it_pipe])\n", - " LA_ist_vec[it_pipe] = T1.get_current_LA()\n", + " turbine.update_LA(LA_soll_vec[it_pipe])\n", + " LA_ist_vec[it_pipe] = turbine.get_current_LA()\n", "\n", - " T1.set_pressure(p_old[-1])\n", - " # set boundary conditions for the next timestep of the characteristic method\n", - " p_boundary_res[it_pipe] = V.get_current_pressure()\n", - " v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_current_Q()\n", + " # set boundary condition for the next timestep of the characterisCon_T_ic method\n", + " turbine.set_pressure(p_old[-1])\n", + " convergence_parameters[0] = p_old[-2]\n", + " convergence_parameters[1] = v_old[-2]\n", + " turbine.converge(convergence_parameters)\n", + " p_boundary_res[it_pipe] = reservoir.get_current_pressure()\n", + " v_boundary_tur[it_pipe] = 1/Pip_area*turbine.get_current_Q()\n", + " Q_boundary_tur[it_pipe] = turbine.get_current_Q()\n", "\n", - " # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n", + " # the the boundary condition in the pipe.object and thereby calculate boundary pressure at turbine\n", " pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n", + " pipe.v[0] = (0.8*pipe.v[0]+0.2*reservoir.get_current_outflux()/Res_area_out)\n", " p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\n", " v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n", + " Q_boundary_res[it_pipe] = pipe.get_current_flux_distribution()[0]\n", "\n", - " error_vec[it_pipe] = abs(v_boundary_res[it_pipe]-V.get_current_outflux()/A_pipe)\n", - "\n", - " # perform the next timestep via the characteristic method\n", + " # perform the next timestep via the characterisCon_T_ic method\n", " pipe.timestep_characteristic_method()\n", "\n", " # prepare for next loop\n", " p_old = pipe.get_current_pressure_distribution()\n", " v_old = pipe.get_current_velocity_distribution()\n", + " Q_old = pipe.get_current_flux_distribution()\n", + "\n", "\n", " # plot some stuff\n", " # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n", - " lo_00.remove()\n", - " lo_01.remove()\n", - " # lo_02.remove()\n", + " lo_p.remove()\n", + " lo_pmin.remove()\n", + " lo_pmax.remove()\n", + " lo_q.remove()\n", + " lo_qmin.remove()\n", + " lo_qmax.remove()\n", " # plot new pressure and velocity distribution in the pipeline\n", - " lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.',c='blue')\n", - " lo_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n", - " # lo_02, = axs1[2].plot(level_vec_2,c='blue')\n", + " lo_p, = axs1[0].plot(Pip_x_vec,pipe.get_current_pressure_distribution(disp=True),marker='.',c='blue')\n", + " lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n", + " lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", + " lo_q, = axs1[1].plot(Pip_x_vec,pipe.get_current_flux_distribution(),marker='.',c='blue')\n", + " lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n", + " lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n", " fig1.suptitle(str(round(t_vec[it_pipe],2))+ ' s / '+str(round(t_vec[-1],2)) + ' s' )\n", " fig1.canvas.draw()\n", " fig1.tight_layout()\n", " fig1.show()\n", - " plt.pause(0.001) \n", - "\n", - " \n", - " " + " plt.pause(0.001) " ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 13, "metadata": {}, "outputs": [], "source": [ - "# plot time evolution of boundary pressure and velocity as well as the reservoir level\n", + "# plot Con_T_ime evolution of boundary pressure and velocity as well as the reservoir level\n", "\n", - "fig2,axs2 = plt.subplots(3,2)\n", - "axs2[0,0].set_title('Pressure reservoir')\n", - "axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,initial_pressure_unit, conversion_pressure_unit))\n", - "axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[0,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", + "fig2,axs2 = plt.subplots(1,1)\n", + "axs2.set_title('Level and Volume reservoir')\n", + "axs2.plot(t_vec,level_vec,label='level')\n", + "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", + "axs2.set_ylabel(r'$h$ [m]')\n", + "x_twin_00 = axs2.twinx()\n", + "x_twin_00.set_ylabel(r'$V$ [$\\mathrm{m}^3$]')\n", + "x_twin_00.plot(t_vec,volume_vec)\n", + "axs2.legend()\n", "\n", - "axs2[0,1].set_title('Velocity reservoir')\n", - "axs2[0,1].plot(t_vec,v_boundary_res)\n", - "axs2[0,1].set_ylim(-2*Q_nenn,+2*Q_nenn)\n", - "axs2[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", + "fig2,axs2 = plt.subplots(1,1)\n", + "axs2.set_title('LA')\n", + "axs2.plot(t_vec,100*LA_soll_vec,label='Target')\n", + "axs2.plot(t_vec,100*LA_ist_vec,label='Actual')\n", + "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", + "axs2.set_ylabel(r'$LA$ [%]')\n", + "axs2.legend()\n", "\n", - "axs2[1,0].set_title('Pressure turbine')\n", - "axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,initial_pressure_unit, conversion_pressure_unit))\n", - "axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[1,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n", + "fig2,axs2 = plt.subplots(1,1)\n", + "axs2.set_title('Pressure reservoir and turbine')\n", + "axs2.plot(t_vec,pressure_conversion(p_boundary_res,pUnit_calc, pUnit_conv),label='Reservoir')\n", + "axs2.plot(t_vec,pressure_conversion(p_boundary_tur,pUnit_calc, pUnit_conv),label='Turbine')\n", + "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", + "axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n", + "axs2.legend()\n", "\n", - "axs2[1,1].set_title('Velocity turbine')\n", - "axs2[1,1].plot(t_vec,v_boundary_tur)\n", - "axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", + "fig2,axs2 = plt.subplots(1,1)\n", + "axs2.set_title('Fluxes')\n", + "axs2.plot(t_vec,Q_boundary_res,label='Outflux')\n", + "axs2.plot(t_vec,Q_in_vec,label='Influx')\n", + "axs2.plot(t_vec,Q_boundary_tur,label='Flux Turbine')\n", + "axs2.set_ylim(-2*Tur_Q_nenn,+2*Tur_Q_nenn)\n", + "axs2.set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", + "axs2.set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n", + "axs2.legend()\n", + "\n", + "fig2,axs2 = plt.subplots(1,1)\n", + "axs2.set_title('Min and Max Pressure')\n", + "axs2.plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n", + "axs2.plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", + "axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", + "axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n", + "\n", + "fig2,axs2 = plt.subplots(1,1)\n", + "axs2.set_title('Min and Max Fluxes')\n", + "axs2.plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n", + "axs2.plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n", + "axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", + "axs2.set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n", + "\n", + "# axs2[0,1].legend()\n", + "# axs2[1,0].legend()\n", + "# axs2[1,1].legend()\n", + "# # axs2[2,0].legend()\n", + "# # axs2[2,1].legend()\n", "\n", - "axs2[2,0].set_title('Level reservoir')\n", - "axs2[2,0].plot(t_vec,level_vec)\n", - "axs2[2,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[2,0].set_ylabel(r'$h$ [m]')\n", "\n", - "axs2[2,1].set_title('LA')\n", - "axs2[2,1].plot(t_vec,100*LA_soll_vec)\n", - "axs2[2,1].plot(t_vec,100*LA_ist_vec)\n", - "axs2[2,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs2[2,1].set_ylabel(r'$LA$ [%]')\n", "fig2.tight_layout()\n", "plt.show()" ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[]" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "plt.semilogy(t_vec,error_vec)" - ] } ], "metadata": { diff --git a/untertweng.txt b/untertweng.txt deleted file mode 100644 index 1c6adb3..0000000 --- a/untertweng.txt +++ /dev/null @@ -1,22 +0,0 @@ -L = 535 m dn 800 mm -478 m dn 1000 mm -Ersatzdurchmesser - -h_pipe - -h 851.78 Pegel + Leitungsgefälle -Leitungsgefälle: 113 - -Fläche 4.25x10.5 + 30m² = 74 m² -Pegelminimum: 851.18 m - -Unterwasserpegel 738.56 -Gesamtfallhöhe = 851.78-738.56 - -Rohrreibung: 0.014 f_D = lambda -c = 500 m/s - -Q_0 = 100%*0.75+30%*0.75 -Q_extrem = 30%*0.75 - -Q = LA*Q_nenn*sqrt(H/H_n)