added Druckrohrleitungs class, based on ETH Code

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)

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",

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 update_volume(self):
-    def set_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

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()        
+
+
+
+
+
+
+
+
+
+
+

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"
    }
   }
  },

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
+}