diff --git a/Ausgleichsbecken/dynamic_pipeline_pressure/Ausgleichsbecken_class_file.py b/Ausgleichsbecken/dynamic_pipeline_pressure/Ausgleichsbecken_class_file.py index 64716ea..750f952 100644 --- a/Ausgleichsbecken/dynamic_pipeline_pressure/Ausgleichsbecken_class_file.py +++ b/Ausgleichsbecken/dynamic_pipeline_pressure/Ausgleichsbecken_class_file.py @@ -4,11 +4,11 @@ class Ausgleichsbecken_class: # units area_unit = r'$\mathrm{m}^2$' area_outflux_unit = r'$\mathrm{m}^2$' - level_unit = 'm' - volume_unit = r'$\mathrm{m}^3$' flux_unit = r'$\mathrm{m}^3/\mathrm{s}$' - time_unit = 's' + level_unit = 'm' pressure_unit = 'Pa' + time_unit = 's' + volume_unit = r'$\mathrm{m}^3$' # init def __init__(self,area,outflux_area,level_min,level_max,timestep = 1): @@ -73,6 +73,7 @@ class Ausgleichsbecken_class: h = self.level dt = self.timestep p,_ = pressure_conversion(self.pressure,self.pressure_unit,'Pa') + # update to include p_halfstep p_hs,_ = pressure_conversion(self.pressure,self.pressure_unit,'Pa') alpha = (self.area_outflux/self.area-1) h_hs = self.update_level(dt/2) diff --git a/Ausgleichsbecken/dynamic_pipeline_pressure/Main_Program.ipynb b/Ausgleichsbecken/dynamic_pipeline_pressure/Main_Program.ipynb index 233f368..d127356 100644 --- a/Ausgleichsbecken/dynamic_pipeline_pressure/Main_Program.ipynb +++ b/Ausgleichsbecken/dynamic_pipeline_pressure/Main_Program.ipynb @@ -91,6 +91,7 @@ "i_max = -1\n", "\n", "for i in range(np.size(time_vec)-1):\n", + " # update to include p_halfstep\n", " V.pressure = pressure_vec[i]\n", " V.e_RK_4()\n", " V.level = V.update_level(V.timestep)\n", diff --git a/Ausgleichsbecken/static_pipeline_pressure/Ausgleichsbecken_class_file.py b/Ausgleichsbecken/static_pipeline_pressure/Ausgleichsbecken_class_file.py index 969dc04..b6ce3d5 100644 --- a/Ausgleichsbecken/static_pipeline_pressure/Ausgleichsbecken_class_file.py +++ b/Ausgleichsbecken/static_pipeline_pressure/Ausgleichsbecken_class_file.py @@ -18,13 +18,13 @@ class Ausgleichsbecken_class: self.level_max = level_max # highest allowed water level self.timestep = timestep # timestep of the simulation -# setter - def set_volume(self): + def update_volume(self): self.volume = self.level*self.area +# setter def set_initial_level(self,initial_level): self.level = initial_level - self.set_volume() + self.update_volume() def set_influx(self,influx): self.influx = influx @@ -61,7 +61,10 @@ class Ausgleichsbecken_class: print('The current outflux is', self.outflux, self.flux_unit) # methods + + def update_level(self,timestep): + # dont update volume here, because update_level gets called to calculate h_halfstep net_flux = self.influx-self.outflux delta_V = net_flux*timestep new_level = (self.volume+delta_V)/self.area diff --git a/Druckrohrleitung/Druckrohrleitung_class_file.py b/Druckrohrleitung/Druckrohrleitung_class_file.py new file mode 100644 index 0000000..8b81e9f --- /dev/null +++ b/Druckrohrleitung/Druckrohrleitung_class_file.py @@ -0,0 +1,162 @@ +from pressure_conversion import pressure_conversion +import numpy as np + +class Druckrohrleitung_class: +# units + acceleration_unit = r'$\mathrm{m}/\mathrm{s}^2$' + angle_unit = '°' + area_unit = r'$\mathrm{m}^2$' + density_unit = r'$\mathrm{kg}/\mathrm{m}^3$' + flux_unit = r'$\mathrm{m}^3/\mathrm{s}$' + length_unit = 'm' + pressure_unit = 'Pa' + time_unit = 's' + velocity_unit = r'$\mathrm{m}/\mathrm{s}$' # for flux and pressure propagation + volume_unit = r'$\mathrm{m}^3$' + + +# init + def __init__(self,total_length,diameter,number_segments,pipeline_angle,Darcy_friction_factor,rho=1000,g=9.81): + self.length = total_length + self.dia = diameter + self.n_seg = number_segments + self.angle = pipeline_angle + self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient + self.density = 1000 + self.g = g + + self.dx = total_length/number_segments + self.l_vec = np.arange(0,(number_segments+1)*self.dx,self.dx) + + # workaround for try-except construct in set_number_of_timesteps + self.c = 0 + +# setter + def set_pressure_propagation_velocity(self,c): + self.c = c + self.dt = self.dx/c + + def set_number_of_timesteps(self,number_timesteps): + self.nt = number_timesteps + if self.c == 0: + 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,input_unit = 'Pa'): + p,_ = pressure_conversion(pressure,input_unit,target_unit=self.pressure_unit) + if np.size(p) == 1: + self.p0 = np.full_like(self.l_vec,p) + elif np.size(p) == np.size(self.l_vec): + self.p0 = p + else: + raise Exception('Unable to assign initial pressure. Input has to be of size 1 or' + np.size(self.l_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_new = np.empty_like(self.p_old) + + def set_initial_flow_velocity(self,velocity): + 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 + else: + raise Exception('Unable to assign initial velocity. Input has to be of size 1 or' + np.size(self.l_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_new = np.empty_like(self.v_old) + + def set_boundary_conditions_next_timestep(self,v_reservoir,p_reservoir,v_turbine,input_unit_pressure = 'Pa'): + rho = self.density + c = self.c + f_D = self.f_D + dt = self.dt + D = self.dia + p_old = self.p_old[-2] # @ second to last node (the one before the turbine) + v_old = self.v_old[-2] # @ second to last node (the one before the turbine) + self.v_boundary_res = v_reservoir + self.v_boundary_tur = v_turbine + self.p_boundary_res,_ = pressure_conversion(p_reservoir,input_unit_pressure,target_unit=self.pressure_unit) + self.p_boundary_tur = p_old+rho*c*v_old-rho*c*f_D*dt/(2*D)*abs(v_old)*v_old + self.v_new[0] = self.v_boundary_res.copy() + self.v_new[-1] = self.v_boundary_tur.copy() + self.p_new[0] = self.p_boundary_res.copy() + self.p_new[-1] = self.p_boundary_tur.copy() + +# getter + def get_pipeline_geometry(self): + print('The total length of the pipeline is', '\n', \ + self.length, self.length_unit, '\n', \ + 'The diameter of the pipeline is', '\n', \ + self.dia, self.length_unit, '\n', \ + 'The pipeline is divided into', self.n_seg , 'segments of length', '\n', \ + round(self.dx,1), self.length_unit, '\n', \ + 'The pipeline has an inclination angle of', '\n', \ + self.angle, self.angle_unit) + + def get_other_pipeline_info(self): + print('The Darcy-friction factor of the pipeline is', '\n', \ + self.f_D, '\n', \ + 'The pipeline is filled with a liquid with density', '\n', \ + self.density, self.density_unit, '\n', \ + 'The gravitational acceleration is set to', '\n', \ + self.g, self.acceleration_unit) + + def get_pressure_propagation_velocity(self): + print('The pressure propagation velocity in the pipeline is', '\n', \ + self.c, self.velocity_unit) + + def get_number_of_timesteps(self): + print(self.nt, 'timesteps are performed in the simulation') + + + def get_initial_pressure(self,target_unit='bar'): + print('The inital pressure distribution in is', '\n', \ + pressure_conversion(self.p0,self.pressure_unit,target_unit)) + + def get_initial_flow_velocity(self): + print('The inital velocity distribution is', '\n', \ + self.v0, self.velocity_unit) + + def get_boundary_conditions_next_timestep(self,target_unit_pressure ='bar'): + print('The pressure at the reservoir for the next timestep is', '\n', \ + pressure_conversion(self.p_boundary_res,self.pressure_unit,target_unit_pressure), '\n', \ + 'The velocity at the reservoir for the next timestep is', '\n', \ + self.v_boundary_res, self.velocity_unit, '\n', \ + 'The pressure at the turbine for the next timestep is', '\n', \ + pressure_conversion(self.p_boundary_tur,self.pressure_unit,target_unit_pressure), '\n', \ + 'The velocity at the turbine for the next timestep is', '\n', \ + self.v_boundary_tur, self.velocity_unit) + + + def timestep_characteristic_method(self): + #number of nodes + nn = self.n_seg+1 + rho = self.density + c = self.c + f_D = self.f_D + dt = self.dt + D = self.dia + + for i in range(1,nn-1): + self.v_new[i] = 0.5*(self.v_old[i-1]+self.v_old[i+1])+0.5/(rho*c)*(self.p_old[i-1]-self.p_old[i+1]) \ + -f_D*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]) + + self.p_new[i] = 0.5*rho*c*(self.v_old[i-1]-self.v_old[i+1])+0.5*(self.p_old[i-1]+self.p_old[i+1]) \ + -rho*c*f_D*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]) + + self.p_old = self.p_new.copy() + self.v_old = self.v_new.copy() + + + + + + + + + + + diff --git a/Druckrohrleitung/Druckrohrleitung.py b/Druckrohrleitung/Druckrohrleitung_functions.py similarity index 100% rename from Druckrohrleitung/Druckrohrleitung.py rename to Druckrohrleitung/Druckrohrleitung_functions.py diff --git a/Druckrohrleitung/Druckstoß_ETH.ipynb b/Druckrohrleitung/Druckstoß_ETH.ipynb index 2251655..fafcd56 100644 --- a/Druckrohrleitung/Druckstoß_ETH.ipynb +++ b/Druckrohrleitung/Druckstoß_ETH.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 22, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -14,7 +14,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -28,14 +28,14 @@ "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 = 1500 # number of time steps after initial conditions\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]" ] }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -70,7 +70,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -90,7 +90,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -139,7 +139,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -159,7 +159,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3.9.7 ('base')", + "display_name": "Python 3.8.13 ('Georg_DT_Slot3')", "language": "python", "name": "python3" }, @@ -173,12 +173,12 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.7" + "version": "3.8.13" }, "orig_nbformat": 4, "vscode": { "interpreter": { - "hash": "ad2bdc8ecc057115af97d19610ffacc2b4e99fae6737bb82f5d7fb13d2f2c186" + "hash": "84fb123bdc47ab647d3782661abcbe80fbb79236dd2f8adf4cef30e8755eb2cd" } } }, diff --git a/Druckrohrleitung/Main_Programm.ipynb b/Druckrohrleitung/Main_Programm.ipynb new file mode 100644 index 0000000..b6deb3f --- /dev/null +++ b/Druckrohrleitung/Main_Programm.ipynb @@ -0,0 +1,218 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from Druckrohrleitung_class_file import Druckrohrleitung_class\n", + "import matplotlib.pyplot as plt\n", + "from pressure_conversion import pressure_conversion" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "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": null, + "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_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": null, + "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", + "\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_new)\n", + " lo_01.set_ydata(pipe.v_new)\n", + "\n", + " # store parameters of node 1 (at reservoir)\n", + " pipe.p_1[it] = pipe.p_new[0]\n", + " pipe.v_1[it] = pipe.v_new[0]\n", + " # store parameters of node N+1 (at reservoir)\n", + " pipe.p_np1[it] = pipe.p_new[-1]\n", + " pipe.v_np1[it] = pipe.v_new[-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": null, + "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", + "fig3.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 +}