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

working on a fix for steady state Ausgleichsbecken

Browse Source
This commit is contained in:
Brantegger Georg
2022-07-27 16:02:39 +02:00
parent 176fa556ff
commit 03ff67e0ad
4 changed files with 121 additions and 46 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Ausgleichsbecken/Ausgleichsbecken_class_file.py
View File

@@ -21,7 +21,7 @@ def FODE_function(x,h,A,A_a,p,rho,g):
# A ... Reservoir_Area # A ... Reservoir_Area
# g ... gravitational acceleration # g ... gravitational acceleration
# rho ... density of the liquid in the reservoir # 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 return f

97
Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 1, "execution_count": 7,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -21,32 +21,64 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 2, "execution_count": 8,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "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", "# define constants\n",
"initial_level = 10. # m\n", "initial_level = 10.1 # m\n",
"initial_influx = 5. # m³/s\n", "initial_influx = 0.8 # m³/s\n",
"# initial_outflux = 1. # m³/s\n",
"# initial_pipeline_pressure = 10.\n",
"# initial_pressure_unit = 'mWS'\n",
"conversion_pressure_unit = 'mWS'\n", "conversion_pressure_unit = 'mWS'\n",
"\n", "\n",
"area_base = 1. # m²\n", "area_base = 75. # m²\n",
"area_outflux = 0.5 # m²\n", "area_outflux = (0.9/2)**2*np.pi # m²\n",
"critical_level_low = 0. # m\n", "critical_level_low = 0. # m\n",
"critical_level_high = 10. # m\n", "critical_level_high = 10. # m\n",
"simulation_timestep = 0.001 # s\n", "simulation_timestep = dt # s\n",
"\n", "\n",
"# for while loop\n", "# for while loop\n",
"total_min_level = 0.01 # m\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", "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": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -60,22 +92,25 @@
"# V.pressure = converted_pressure\n", "# V.pressure = converted_pressure\n",
"V.set_steady_state(initial_influx,initial_level,conversion_pressure_unit)\n", "V.set_steady_state(initial_influx,initial_level,conversion_pressure_unit)\n",
"\n", "\n",
"time_vec = np.arange(0,total_max_time,simulation_timestep)\n", "time_vec = np.arange(0,nt+1,1)*simulation_timestep\n",
"outflux_vec = np.empty_like(time_vec)\n", "outflux_vec = np.zeros_like(time_vec)\n",
"outflux_vec[0] = V.get_current_outflux()\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", "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", "\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,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", " \n",
"i_max = -1\n", "i_max = -1\n",
"\n", "\n",
"for i in range(np.size(time_vec)-1):\n", "for i in range(1,nt+1):\n",
" V.set_pressure(pressure_vec[i])\n", " V.set_pressure(pressure_vec[i-1])\n",
" V.set_outflux(outflux_vec[i-1])\n",
" V.timestep_reservoir_evolution()\n", " V.timestep_reservoir_evolution()\n",
" outflux_vec[i+1] = V.get_current_outflux()\n", " outflux_vec[i] = V.get_current_outflux()\n",
" level_vec[i+1] = V.get_current_level()\n", " level_vec[i] = V.get_current_level()\n",
" pressure_vec[i] = V.get_current_pressure()\n",
" if V.level < total_min_level:\n", " if V.level < total_min_level:\n",
" i_max = i\n", " i_max = i\n",
" break\n", " break\n",
@@ -84,7 +119,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 7, "execution_count": 11,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -111,6 +146,26 @@
"\n", "\n",
"fig1.tight_layout() " "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": { "metadata": {

2
Druckrohrleitung/Druckrohrleitung_class_file.py
View File

@@ -117,7 +117,7 @@ class Druckrohrleitung_class:
# the flow velocity is given by the constant flow through the pipe # 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(self.n_seg+1,ss_flux/self.A)
# the static pressure is given by the hydrostatic pressure, corrected for friction losses and dynamic pressure # 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_flow_velocity(ss_v0)
self.set_initial_pressure(ss_pressure) self.set_initial_pressure(ss_pressure)

44
Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb
View File

@@ -35,8 +35,8 @@
"L = 1000. # length of pipeline [m]\n", "L = 1000. # length of pipeline [m]\n",
"D = 0.9 # pipe diameter [m]\n", "D = 0.9 # pipe diameter [m]\n",
"h_res = 10. # water level in upstream reservoir [m]\n", "h_res = 10. # water level in upstream reservoir [m]\n",
"n = 50 # number of pipe segments in discretization\n", "n = 50000 # number of pipe segments in discretization\n",
"nt = 5000 # number of time steps after initial conditions\n", "nt = 12 # number of time steps after initial conditions\n",
"f_D = 0.01 # Darcy friction factor\n", "f_D = 0.01 # Darcy friction factor\n",
"c = 400. # propagation velocity of the pressure wave [m/s]\n", "c = 400. # propagation velocity of the pressure wave [m/s]\n",
"h_pipe = 105. # hydraulic head without reservoir [m] \n", "h_pipe = 105. # hydraulic head without reservoir [m] \n",
@@ -49,9 +49,9 @@
"dx = L/n # length of each pipe segment\n", "dx = L/n # length of each pipe segment\n",
"dt = dx/c # timestep according to method of characterisitics\n", "dt = dx/c # timestep according to method of characterisitics\n",
"nn = n+1 # number of nodes\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", "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*dt,dt) # time vector\n", "t_vec = np.arange(0,nt,1)*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", "h_vec = np.arange(0,nn,1)*h_pipe/n # hydraulic head of pipeline at each node\n",
"\n", "\n",
"\n", "\n",
"# define constants reservoir\n", "# define constants reservoir\n",
@@ -91,6 +91,7 @@
"print(V.get_current_influx())\n", "print(V.get_current_influx())\n",
"print(V.get_current_outflux())\n", "print(V.get_current_outflux())\n",
"print(V.get_current_level())\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(V.get_current_pressure())\n",
"print(pipe.get_current_pressure_distribution()[0])\n", "print(pipe.get_current_pressure_distribution()[0])\n",
"print(pipe.get_current_velocity_distribution()*area_pipe)\n", "print(pipe.get_current_velocity_distribution()*area_pipe)\n",
@@ -99,7 +100,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 5, "execution_count": null,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -130,7 +131,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 6, "execution_count": null,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -153,16 +154,33 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 7, "execution_count": 22,
"metadata": {}, "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<cell line: 1>\u001b[1;34m()\u001b[0m\n\u001b[0;32m <a href='vscode-notebook-cell:/y%3A/KELAG/KS/KS-PW/04%20Digitalisierung/KSPWDEV%20Server/Digital%20Trainee%20Projekt/DT_Slot_3_Project_Repo/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb#ch0000006?line=39'>40</a>\u001b[0m fig1\u001b[39m.\u001b[39mcanvas\u001b[39m.\u001b[39mdraw()\n\u001b[0;32m <a href='vscode-notebook-cell:/y%3A/KELAG/KS/KS-PW/04%20Digitalisierung/KSPWDEV%20Server/Digital%20Trainee%20Projekt/DT_Slot_3_Project_Repo/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb#ch0000006?line=40'>41</a>\u001b[0m fig1\u001b[39m.\u001b[39mtight_layout()\n\u001b[1;32m---> <a href='vscode-notebook-cell:/y%3A/KELAG/KS/KS-PW/04%20Digitalisierung/KSPWDEV%20Server/Digital%20Trainee%20Projekt/DT_Slot_3_Project_Repo/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb#ch0000006?line=41'>42</a>\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": [ "source": [
"\n", "\n",
"for it_pipe in range(1,nt):\n", "for it_pipe in range(1,nt):\n",
"# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\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", " # set initial conditions for the reservoir time evolution calculted with e-RK4\n",
" V.set_pressure = p_old[0]\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", " # 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", " for it_res in range(nt_eRK4):\n",
" V.timestep_reservoir_evolution() \n", " V.timestep_reservoir_evolution() \n",
@@ -171,6 +189,7 @@
" \n", " \n",
" # set boundary conditions for the next timestep of the characteristic method\n", " # set boundary conditions for the next timestep of the characteristic method\n",
" p_boundary_res[it_pipe] = V.get_current_pressure()\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", " v_boundary_tur[it_pipe] = initial_flux/area_pipe\n",
"\n", "\n",
" # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\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_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_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n",
" \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.canvas.draw()\n",
" fig1.tight_layout()\n", " fig1.tight_layout()\n",
" plt.pause(0.00001)\n" " plt.pause(0.000001)\n",
"\n"
] ]
}, },
{ {