{ "cells": [ { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from Druckrohrleitung_class_file import Druckrohrleitung_class\n", "import matplotlib.pyplot as plt\n", "\n", "#importing pressure conversion function\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" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "%matplotlib qt5\n", "#define constants\n", "\n", "g = 9.81 # gravitational acceleration [m/s²]\n", "\n", "L = 1000 # length of pipeline [m]\n", "rho = 1000 # density of water [kg/m³]\n", "D = 1 # pipe diameter [m]\n", "Q0 = 2 # initial flow in whole pipe [m³/s]\n", "h = 20 # water level in upstream reservoir [m]\n", "n = 10 # number of pipe segments in discretization\n", "nt = 500 # 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", "\n", "\n", "# preparing the discretization and initial conditions\n", "\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*dx,dx) # pl = pipe-length. position of the nodes on the pipeline\n", "t_vec = np.arange(0,nt*dt,dt) # time vector\n", "\n", "v0 = Q0/(D**2/4*np.pi)\n", "p0 = (rho*g*h-v0**2*rho/2)\n", "\n", "# storage vectors for old parameters\n", "v_old = np.full(nn,v0)\n", "p_old = p0-(f_D*pl_vec/D*rho/2*v0**2) # ref Wikipedia: Darcy Weisbach\n", "\n", "# storage vectors for new parameters\n", "v_new = np.zeros_like(v_old)\n", "p_new = np.zeros_like(p_old)\n", "\n", "# storage vector for time evolution of parameters at node 1 (at reservoir)\n", "p_1 = np.full_like(t_vec,p0)\n", "v_1 = np.full_like(t_vec,v0)\n", "\n", "# storage vector for time evolution of parameters at node N+1 (at valve)\n", "p_np1 = np.full_like(t_vec,p0)\n", "v_np1 = np.full_like(t_vec,v0)\n", "\n", "for it in range(1,nt):\n", "\n", " # set boundary conditions\n", " v_new[-1] = 0 # in front of the instantaneously closing valve, the velocity is 0\n", " p_new[0] = p0 # hydrostatic pressure from the reservoir\n", "\n", " # calculate the new parameters at first and last node\n", " v_new[0] = v_old[1]+1/(rho*c)*(p0-p_old[1])-f_D*dt/(2*D)*abs(v_old[1])*v_old[1]\n", " p_new[-1] = p_old[-2]+rho*c*v_old[-2]-rho*c*f_D*dt/(2*D) *abs(v_old[-2])*v_old[-2]\n", "\n", " # calculate parameters at second to second-to-last nodes \n", " #equation 2-30 plus 2-31 (and refactor for v_i^j+1) in block 2\n", "\n", " for i in range(1,nn-1):\n", " v_new[i] = 0.5*(v_old[i-1]+v_old[i+1])+0.5/(rho*c)*(p_old[i-1]-p_old[i+1]) \\\n", " -f_D*dt/(4*D)*(abs(v_old[i-1])*v_old[i-1]+abs(v_old[i+1])*v_old[i+1])\n", "\n", " p_new[i] = 0.5*rho*c*(v_old[i-1]-v_old[i+1])+0.5*(p_old[i-1]+p_old[i+1]) \\\n", " -rho*c*f_D*dt/(4*D)*(abs(v_old[i-1])*v_old[i-1]-abs(v_old[i+1])*v_old[i+1])\n", " \n", "\n", " # prepare for next loop\n", " # use .copy() to avoid that memory address is overwritten and hell breaks loose :D\n", " #https://www.geeksforgeeks.org/array-copying-in-python/\n", " p_old = p_new.copy()\n", " v_old = v_new.copy()\n", "\n", " # store parameters of node 1 (at reservoir)\n", " p_1[it] = p_new[0]\n", " v_1[it] = v_new[0]\n", " # store parameters of node N+1 (at reservoir)\n", " p_np1[it] = p_new[-1]\n", " v_np1[it] = v_new[-1]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "fig1,axs1 = plt.subplots(2,2)\n", "axs1[0,0].plot(t_vec,p_1)\n", "axs1[0,1].plot(t_vec,v_1)\n", "axs1[1,0].plot(t_vec,p_np1)\n", "axs1[1,1].plot(t_vec,v_np1)\n", "axs1[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs1[0,0].set_ylabel(r'$p$ [Pa]')\n", "axs1[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs1[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", "axs1[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs1[1,0].set_ylabel(r'$p$ [Pa]')\n", "axs1[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs1[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", "\n", "axs1[0,0].set_title('Pressure Reservoir')\n", "axs1[0,1].set_title('Velocity Reservoir')\n", "axs1[1,0].set_title('Pressure Turbine')\n", "axs1[1,1].set_title('Velocity Turbine')\n", "fig1.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "pipe = Druckrohrleitung_class(L,D,n,0,f_D)\n", "\n", "pipe.set_pressure_propagation_velocity(c)\n", "pipe.set_number_of_timesteps(nt)\n", "\n", "pipe.set_initial_pressure(p0)\n", "pipe.set_initial_flow_velocity(v0)\n", "pipe.set_boundary_conditions_next_timestep(v_1[0],p_1[0],v_np1[0])\n", "\n", "# storage vector for time evolution of parameters at node 1 (at reservoir)\n", "pipe.p_1 = np.full_like(t_vec,p0)\n", "pipe.v_1 = np.full_like(t_vec,v0)\n", "\n", "# storage vector for time evolution of parameters at node N+1 (at valve)\n", "pipe.p_np1 = np.full_like(t_vec,p0)\n", "pipe.v_np1 = np.full_like(t_vec,v0)\n", "\n", "fig2,axs2 = plt.subplots(2,1)\n", "axs2[0].set_title('Pressure distribution in pipeline')\n", "axs2[1].set_title('Velocity distribution in pipeline')\n", "axs2[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs2[0].set_ylabel(r'$p$ [Pa]')\n", "axs2[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs2[1].set_ylabel(r'$p$ [Pa]')\n", "lo_00, = axs2[0].plot(pl_vec,pipe.p_old,marker='.')\n", "lo_01, = axs2[1].plot(pl_vec,pipe.v_old,marker='.')\n", "axs2[0].set_ylim([-20*p0,20*p0])\n", "axs2[1].set_ylim([-2*v0,2*v0])\n", "fig2.tight_layout()\n", "\n", "\n", "for it in range(1,pipe.nt):\n", " pipe.set_boundary_conditions_next_timestep(v_1[it],p_1[it],v_np1[it])\n", " pipe.timestep_characteristic_method()\n", " lo_00.set_ydata(pipe.p)\n", " lo_01.set_ydata(pipe.v)\n", "\n", " # store parameters of node 1 (at reservoir)\n", " pipe.p_1[it] = pipe.p[0]\n", " pipe.v_1[it] = pipe.v[0]\n", " # store parameters of node N+1 (at reservoir)\n", " pipe.p_np1[it] = pipe.p[-1]\n", " pipe.v_np1[it] = pipe.v[-1]\n", " \n", " fig2.suptitle(str(it))\n", " fig2.canvas.draw()\n", " fig2.tight_layout()\n", " plt.pause(0.001)\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "fig3,axs3 = plt.subplots(2,2)\n", "axs3[0,0].plot(t_vec,pipe.p_1)\n", "axs3[0,1].plot(t_vec,pipe.v_1)\n", "axs3[1,0].plot(t_vec,pipe.p_np1)\n", "axs3[1,1].plot(t_vec,pipe.v_np1)\n", "axs3[0,0].set_title('Pressure Reservoir')\n", "axs3[0,1].set_title('Velocity Reservoir')\n", "axs3[1,0].set_title('Pressure Turbine')\n", "axs3[1,1].set_title('Velocity Turbine')\n", "axs3[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs3[0,0].set_ylabel(r'$p$ [Pa]')\n", "axs3[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs3[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", "axs3[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs3[1,0].set_ylabel(r'$p$ [Pa]')\n", "axs3[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs3[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", "fig3.tight_layout()\n", "plt.show()" ] } ], "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 }