From 03ff67e0ad8baf586e694bab1fc65823fed25f30 Mon Sep 17 00:00:00 2001 From: Brantegger Georg Date: Wed, 27 Jul 2022 16:02:39 +0200 Subject: [PATCH] working on a fix for steady state Ausgleichsbecken --- .../Ausgleichsbecken_class_file.py | 2 +- .../Ausgleichsbecken_test_steady_state.ipynb | 97 +++++++++++++++---- .../Druckrohrleitung_class_file.py | 24 ++--- .../Druckrohrleitung_test_steady_state.ipynb | 44 ++++++--- 4 files changed, 121 insertions(+), 46 deletions(-) diff --git a/Ausgleichsbecken/Ausgleichsbecken_class_file.py b/Ausgleichsbecken/Ausgleichsbecken_class_file.py index 18350e9..e7f7264 100644 --- a/Ausgleichsbecken/Ausgleichsbecken_class_file.py +++ b/Ausgleichsbecken/Ausgleichsbecken_class_file.py @@ -21,7 +21,7 @@ def FODE_function(x,h,A,A_a,p,rho,g): # A ... Reservoir_Area # g ... gravitational acceleration # rho ... density of the liquid in the reservoir - f = x*abs(x)/h*(A_a/A-1.5)+g-p/(rho*h) + f = x*abs(x)/h*(A_a/A-1.)+g-p/(rho*h) return f diff --git a/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb b/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb index 1b82115..6d6c071 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": 1, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -21,32 +21,64 @@ }, { "cell_type": "code", - "execution_count": 2, + "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. # m\n", - "initial_influx = 5. # m³/s\n", - "# initial_outflux = 1. # m³/s\n", - "# initial_pipeline_pressure = 10.\n", - "# initial_pressure_unit = 'mWS'\n", + "initial_level = 10.1 # m\n", + "initial_influx = 0.8 # m³/s\n", "conversion_pressure_unit = 'mWS'\n", "\n", - "area_base = 1. # m²\n", - "area_outflux = 0.5 # m²\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 = 0.001 # s\n", + "simulation_timestep = dt # s\n", "\n", "# for while loop\n", "total_min_level = 0.01 # m\n", - "total_max_time = 1000 # s" + "total_max_time = 1000 # s\n", + "\n", + "nt = int(total_max_time//simulation_timestep)" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "# # define constants\n", + "# initial_level = 10.1 # m\n", + "# initial_influx = 5. # 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 = 1. # m²\n", + "# area_outflux = 0.5 # m²\n", + "# critical_level_low = 0. # m\n", + "# critical_level_high = 10. # m\n", + "# simulation_timestep = 0.0005 # s\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)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -60,22 +92,25 @@ "# V.pressure = converted_pressure\n", "V.set_steady_state(initial_influx,initial_level,conversion_pressure_unit)\n", "\n", - "time_vec = np.arange(0,total_max_time,simulation_timestep)\n", - "outflux_vec = np.empty_like(time_vec)\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.empty_like(time_vec)\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", - "pressure_vec = np.full_like(time_vec,V.get_current_pressure())\n", " \n", "i_max = -1\n", "\n", - "for i in range(np.size(time_vec)-1):\n", - " V.set_pressure(pressure_vec[i])\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+1] = V.get_current_outflux()\n", - " level_vec[i+1] = V.get_current_level()\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", @@ -84,7 +119,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ @@ -111,6 +146,26 @@ "\n", "fig1.tight_layout() " ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "19.987523898552976" + ] + }, + "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 71d553b..7cf4c12 100644 --- a/Druckrohrleitung/Druckrohrleitung_class_file.py +++ b/Druckrohrleitung/Druckrohrleitung_class_file.py @@ -88,17 +88,17 @@ class Druckrohrleitung_class: # the velocity at the reservoir will be calculated using the backward characteristic # constants for a cleaner formula - rho = self.density - c = self.c - f_D = self.f_D - dt = self.dt - D = self.dia - g = self.g - alpha = self.angle - p_old_tur = self.p_old[-2] # @ second to last node (the one before the turbine) - v_old_tur = self.v_old[-2] # @ second to last node (the one before the turbine) - p_old_res = self.p_old[1] # @ second node (the one after the reservoir) - v_old_res = self.v_old[1] # @ second node (the one after the reservoir) + rho = self.density + c = self.c + f_D = self.f_D + dt = self.dt + D = self.dia + g = self.g + alpha = self.angle + p_old_tur = self.p_old[-2] # @ second to last node (the one before the turbine) + v_old_tur = self.v_old[-2] # @ second to last node (the one before the turbine) + p_old_res = self.p_old[1] # @ second node (the one after the reservoir) + v_old_res = self.v_old[1] # @ second node (the one after the reservoir) # set the boundary conditions derived from reservoir and turbine v_boundary_tur = v_turbine # at new timestep p_boundary_res = p_reservoir # at new timestep @@ -117,7 +117,7 @@ class Druckrohrleitung_class: # the flow velocity is given by the constant flow through the pipe ss_v0 = np.full(self.n_seg+1,ss_flux/self.A) # the static pressure is given by the hydrostatic pressure, corrected for friction losses and dynamic pressure - ss_pressure = (self.density*self.g*(ss_level_reservoir+h_vec)-ss_v0**2*self.density/2)-(self.f_D*pl_vec/self.dia*self.density/2*ss_v0**2) + ss_pressure = self.density*self.g*(ss_level_reservoir+h_vec)-ss_v0**2*self.density/2-(self.f_D*pl_vec/self.dia*self.density/2*ss_v0**2) self.set_initial_flow_velocity(ss_v0) self.set_initial_pressure(ss_pressure) diff --git a/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb b/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb index e92c8b6..ff4719f 100644 --- a/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb +++ b/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb @@ -35,8 +35,8 @@ "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", + "n = 50000 # number of pipe segments in discretization\n", + "nt = 12 # 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", @@ -49,9 +49,9 @@ "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*dx,dx) # pl = pipe-length. position of the nodes on the pipeline\n", - "t_vec = np.arange(0,nt*dt,dt) # time vector\n", - "h_vec = np.arange(0,h_pipe+h_pipe/n,h_pipe/n) # hydraulic head of pipeline at each node\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", "# define constants reservoir\n", @@ -91,6 +91,7 @@ "print(V.get_current_influx())\n", "print(V.get_current_outflux())\n", "print(V.get_current_level())\n", + "print(rho*g*V.get_current_level()-rho/2*(V.get_current_outflux()/area_pipe)**2)\n", "print(V.get_current_pressure())\n", "print(pipe.get_current_pressure_distribution()[0])\n", "print(pipe.get_current_velocity_distribution()*area_pipe)\n", @@ -99,7 +100,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -130,7 +131,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -153,16 +154,33 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 22, "metadata": {}, - "outputs": [], + "outputs": [ + { + "ename": "KeyboardInterrupt", + "evalue": "", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", + "\u001b[1;32my:\\KELAG\\KS\\KS-PW\\04 Digitalisierung\\KSPWDEV Server\\Digital Trainee Projekt\\DT_Slot_3_Project_Repo\\Druckrohrleitung\\Druckrohrleitung_test_steady_state.ipynb Cell 7\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 40\u001b[0m fig1\u001b[39m.\u001b[39mcanvas\u001b[39m.\u001b[39mdraw()\n\u001b[0;32m 41\u001b[0m fig1\u001b[39m.\u001b[39mtight_layout()\n\u001b[1;32m---> 42\u001b[0m plt\u001b[39m.\u001b[39;49mpause(\u001b[39m0.000001\u001b[39;49m)\n", + "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\pyplot.py:548\u001b[0m, in \u001b[0;36mpause\u001b[1;34m(interval)\u001b[0m\n\u001b[0;32m 546\u001b[0m canvas\u001b[39m.\u001b[39mdraw_idle()\n\u001b[0;32m 547\u001b[0m show(block\u001b[39m=\u001b[39m\u001b[39mFalse\u001b[39;00m)\n\u001b[1;32m--> 548\u001b[0m canvas\u001b[39m.\u001b[39;49mstart_event_loop(interval)\n\u001b[0;32m 549\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[0;32m 550\u001b[0m time\u001b[39m.\u001b[39msleep(interval)\n", + "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\backends\\backend_qt.py:409\u001b[0m, in \u001b[0;36mFigureCanvasQT.start_event_loop\u001b[1;34m(self, timeout)\u001b[0m\n\u001b[0;32m 405\u001b[0m timer \u001b[39m=\u001b[39m QtCore\u001b[39m.\u001b[39mQTimer\u001b[39m.\u001b[39msingleShot(\u001b[39mint\u001b[39m(timeout \u001b[39m*\u001b[39m \u001b[39m1000\u001b[39m),\n\u001b[0;32m 406\u001b[0m event_loop\u001b[39m.\u001b[39mquit)\n\u001b[0;32m 408\u001b[0m \u001b[39mwith\u001b[39;00m _maybe_allow_interrupt(event_loop):\n\u001b[1;32m--> 409\u001b[0m qt_compat\u001b[39m.\u001b[39m_exec(event_loop)\n", + "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\contextlib.py:120\u001b[0m, in \u001b[0;36m_GeneratorContextManager.__exit__\u001b[1;34m(self, type, value, traceback)\u001b[0m\n\u001b[0;32m 118\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mtype\u001b[39m \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[0;32m 119\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[1;32m--> 120\u001b[0m \u001b[39mnext\u001b[39;49m(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mgen)\n\u001b[0;32m 121\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mStopIteration\u001b[39;00m:\n\u001b[0;32m 122\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mFalse\u001b[39;00m\n", + "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\backends\\qt_compat.py:262\u001b[0m, in \u001b[0;36m_maybe_allow_interrupt\u001b[1;34m(qapp)\u001b[0m\n\u001b[0;32m 260\u001b[0m signal\u001b[39m.\u001b[39msignal(signal\u001b[39m.\u001b[39mSIGINT, old_sigint_handler)\n\u001b[0;32m 261\u001b[0m \u001b[39mif\u001b[39;00m handler_args \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[1;32m--> 262\u001b[0m old_sigint_handler(\u001b[39m*\u001b[39;49mhandler_args)\n", + "\u001b[1;31mKeyboardInterrupt\u001b[0m: " + ] + } + ], "source": [ "\n", "for it_pipe in range(1,nt):\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", + " # V.set_outflux = v_old[0]*area_pipe\n", + " print(V.get_current_pressure())\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", @@ -171,6 +189,7 @@ " \n", " # set boundary conditions for the next timestep of the characteristic method\n", " p_boundary_res[it_pipe] = V.get_current_pressure()\n", + " print(V.get_current_pressure())\n", " v_boundary_tur[it_pipe] = initial_flux/area_pipe\n", "\n", " # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n", @@ -194,10 +213,11 @@ " 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", " \n", - " fig1.suptitle(str(round(t_vec[it_pipe],2)) + '/' + str(t_vec[-1]))\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.00001)\n" + " plt.pause(0.000001)\n", + "\n" ] }, {