Merge branch 'Dev' of https://github.com/GeorgBrantegger/Kelag_DT_Slot_3 into Dev

.gitignore

@@ -1,6 +1,4 @@
 2bignored.txt
-functions/__pycache__/
+*__pycache__/
 .vscode/settings.json
-__pycache__/Ausgleichsbecken_class_file.cpython-38.pyc
-__pycache__/Ausgleichsbecken.cpython-38.pyc
-__pycache__/functions.cpython-38.pyc
+*.pyc

Ausgleichsbecken/dynamic_pipeline_pressure/Ausgleichsbecken_class_file.py

@@ -0,0 +1,87 @@
+from Ausgleichsbecken import FODE_function, get_h_halfstep, get_p_halfstep
+from pressure_conversion import pressure_conversion
+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'
+    pressure_unit       = 'Pa'
+
+# init
+    def __init__(self,area,outflux_area,level_min,level_max,timestep = 1):
+        self.area           = area              # base area of the rectangular structure
+        self.area_outflux   = outflux_area      # area of the outlet towards the pipeline        
+        self.level_min      = level_min         # lowest  allowed water level    
+        self.level_max      = level_max         # highest allowed water level
+        self.timestep       = timestep          # timestep of the simulation    
+
+# setter
+    def set_volume(self):
+        self.volume = self.level*self.area 
+
+    def set_initial_level(self,initial_level):
+        self.level = initial_level
+        self.set_volume()
+
+    def set_influx(self,influx):
+        self.influx = influx
+
+    def set_outflux(self,outflux):
+        self.outflux = outflux
+
+# getter
+    def get_area(self):
+        print('The base area of the cuboid reservoir is', self.area, self.area_unit)
+
+    def get_outflux_area(self):
+        print('The outflux area from the cuboid reservoir to the pipeline is', \
+            self.area_outflux, self.area_outflux_unit)
+    
+    def get_level(self):
+        print('The current level in the reservoir is', self.level , self.level_unit)
+
+    def get_crit_levels(self):
+        print('The critical water levels in the reservoir are:  \n',\
+            '   Minimum:', self.level_min , self.level_unit , '\n',\
+            '   Maximum:', self.level_max , self.level_unit )
+
+    def get_volume(self):
+        print('The current water volume in the reservoir is', self.volume, self.volume_unit)
+
+    def get_timestep(self):
+        print('The timestep for the simulation is' , self.timestep, self.time_unit)
+
+    def get_influx(self):
+        print('The current influx is', self.influx, self.flux_unit)
+
+    def get_outflux(self):
+        print('The current outflux is', self.outflux, self.flux_unit)
+
+# methods
+    def update_level(self,timestep):
+        net_flux = self.influx-self.outflux
+        delta_V = net_flux*timestep
+        new_level = (self.volume+delta_V)/self.area
+        return new_level
+
+
+    def e_RK_4(self):
+        # Update to deal with non constant pipeline pressure!
+        yn = self.outflux/self.area_outflux
+        h = self.level
+        dt = self.timestep
+        p,_ = pressure_conversion(self.pressure,self.pressure_unit,'Pa')
+        p_hs,_ = pressure_conversion(self.pressure,self.pressure_unit,'Pa')
+        alpha = (self.area_outflux/self.area-1)
+        h_hs = self.update_level(dt/2)
+        Y1 = yn
+        Y2 = yn + dt/2*FODE_function(Y1, h, alpha, self.pressure)
+        Y3 = yn + dt/2*FODE_function(Y2, h_hs, alpha, p_hs)
+        Y4 = yn + dt*FODE_function(Y3, h_hs, alpha, p_hs)
+        ynp1 = yn + dt/6*(FODE_function(Y1, h, alpha, p)+2*FODE_function(Y2, h_hs, alpha, p_hs)+ \
+            2*FODE_function(Y3, h_hs, alpha, p_hs)+ FODE_function(Y4, h, alpha, p))
+
+        self.outflux = ynp1*self.area_outflux

Ausgleichsbecken/dynamic_pipeline_pressure/pressure_conversion.py

@@ -0,0 +1,77 @@
+# convert to Pa
+def bar_to_pa(p):
+    return p*1e5
+
+def mWS_to_pa(p):
+    return p*9.80665*1e3
+
+def torr_to_pa(p):
+    return p*133.322
+
+def atm_to_pa(p):
+    return p*101.325*1e3
+
+def psi_to_pa(p):
+    return p*6894.8
+
+# convert from Pa
+def pa_to_bar(p):
+    return p*1e-5
+
+def pa_to_mWS(p):
+    return p*1/(9.80665*1e3)
+
+def pa_to_torr(p):
+    return p/133.322
+
+def pa_to_atm(p):
+    return p*1/(101.325*1e3)
+
+ # converstion function
+
+def pa_to_psi(p):
+    return p/6894.8
+
+def pressure_conversion(pressure, input_unit = 'bar', target_unit = 'Pa'):
+    p = pressure
+    if input_unit.lower()   == 'bar':
+        p_pa = bar_to_pa(p)
+    elif input_unit.lower() == 'mws':
+        p_pa = mWS_to_pa(p)
+    elif input_unit.lower() == 'torr':
+        p_pa = torr_to_pa(p)
+    elif input_unit.lower() == 'atm':
+        p_pa = atm_to_pa(p)
+    elif input_unit.lower() == 'psi':
+        p_pa = psi_to_pa(p)
+    elif input_unit.lower() == 'pa':
+        p_pa = p
+    else:
+        raise Exception('Given input unit not recognised. \n Known units are: Pa, bar, mWs, Torr, atm, psi')
+    
+    if target_unit.lower()      == 'bar':
+        return pa_to_bar(p_pa), target_unit
+    elif target_unit.lower()    == 'mws':
+        return pa_to_mWS(p_pa), target_unit
+    elif target_unit.lower()    == 'torr':
+        return pa_to_torr(p_pa), target_unit
+    elif target_unit.lower()    == 'atm':
+        return pa_to_atm(p_pa), target_unit
+    elif target_unit.lower()    =='psi':
+        return pa_to_psi(p_pa), target_unit
+    elif target_unit.lower()    == 'pa':
+        return p_pa, target_unit
+    else:
+        raise Exception('Given target unit not recognised. \n Known units are: Pa, bar, mWs, Torr, atm, psi')
+
+# testing_pressure_conversion
+if __name__ == '__main__':
+    p = 1
+
+    unit_dict = ['Pa','Bar','Torr','Atm','MWS','psi']
+
+    for input_unit in unit_dict:
+        for target_unit in unit_dict:
+            converted_p = pressure_conversion(p,input_unit,target_unit)
+            print(input_unit,target_unit)
+            print(converted_p)

Ausgleichsbecken/static_pipeline_pressure/Ausgleichsbecken.py

@@ -0,0 +1,67 @@
+import numpy as np 
+
+def Volume_trend(influx, outflux, timestep=1, V_0=0):
+    '''
+    Returns the trend and the volume and the final volume,  defined
+    by influx and outflux patterns. The optional parameter timestep
+    defines the time increment over which the fluxes are changing.
+    '''
+    net_flux = influx-outflux
+    delta_V = net_flux*timestep
+    V_trend = V_0+np.cumsum(delta_V)
+    V_end = V_trend[-1]
+    return V_end,  V_trend
+
+def Height_trend(V_trend, area=1, h_crit_low=-np.inf, h_crit_high=np.inf):
+    '''
+    Returns the trend and the height and the final height,  defined
+    by influx and outflux patterns as well as the crosssection area.
+    The optional parameters h_crit_low/high indicate limits that the height
+    should never exceed. If this occures,  TRUE is returned in the corresponding
+    h_crit_flag.
+    '''
+    h_trend = V_trend/area
+    h_crit_flag_low = np.any(h_trend <= h_crit_low)
+    h_crit_flag_high = np.any(h_trend >= h_crit_high)
+    h_end = h_trend[-1]
+    return h_trend, h_end, h_crit_flag_low, h_crit_flag_high
+
+def get_h_halfstep(initial_height, influx, outflux, timestep, area):
+    h0      = initial_height
+    Q_in    = influx
+    Q_out   = outflux
+    dt      = timestep
+    A       = area
+
+    h_halfstep = h0+1/A*(Q_in-Q_out)*dt/2
+
+def get_p_halfstep(p0, p1):
+    p_halfstep = (p0+p1)/2
+
+def FODE_function(x, h, alpha, p, rho=1000., g=9.81):
+    f = x*abs(x)/h*alpha+g-p/(rho*h)
+    return f
+
+
+def e_RK_4(yn, h, dt, Q0, Q1, A0, A1, p0, p1):
+    alpha = (A1/A0-1)
+
+    h_hs = get_h_halfstep(h, Q0, Q1, dt, A0)
+    p_hs = get_p_halfstep(p0, p1)
+    Y1 = yn
+    Y2 = yn + dt/2*FODE_function(Y1, h, alpha, p0)
+    Y3 = yn + dt/2*FODE_function(Y2, h_hs, alpha, p_hs)
+    Y4 = yn + dt*FODE_function(Y3, h_hs, alpha, p_hs)
+    ynp1 = yn + dt/6*(FODE_function(Y1, h, alpha, p)+2*FODE_function(Y2, h_hs, alpha, p_hs)+ \
+        2*FODE_function(Y3, h_hs, alpha, p_hs)+ FODE_function(Y4, h, alpha, p))
+
+
+
+
+## testing
+# if __name__ == "__main__":
+#     influx = np.full([1, 100],  6)
+#     outflux = np.full_like(influx,  4)
+#     V_end,  V_trend = Volume_trend(influx,  outflux, timestep=0.5, V_0 = 100)
+#     print(V_end)
+#     print(V_trend)

Ausgleichsbecken/static_pipeline_pressure/Ausgleichsbecken_class_file.py

@@ -1,5 +1,5 @@
 from Ausgleichsbecken import FODE_function, get_h_halfstep, get_p_halfstep
-from functions.pressure_conversion import pressure_conversion
+from pressure_conversion import pressure_conversion
 class Ausgleichsbecken_class:
 # units
     area_unit           = r'$\mathrm{m}^2$'

Ausgleichsbecken/static_pipeline_pressure/Main_Program.ipynb

@@ -2,19 +2,19 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
     "import numpy as np\n",
     "from Ausgleichsbecken_class_file import Ausgleichsbecken_class\n",
     "import matplotlib.pyplot as plt\n",
-    "from functions.pressure_conversion import pressure_conversion"
+    "from pressure_conversion import pressure_conversion"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -39,13 +39,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "c7ff9cd7c3ae4d21a32d833650fa73a3",
+       "model_id": "6a4020b97d834285b3b362c8b1f27e47",
        "version_major": 2,
        "version_minor": 0
       },

Ausgleichsbecken/static_pipeline_pressure/pressure_conversion.py

@@ -0,0 +1,77 @@
+# convert to Pa
+def bar_to_pa(p):
+    return p*1e5
+
+def mWS_to_pa(p):
+    return p*9.80665*1e3
+
+def torr_to_pa(p):
+    return p*133.322
+
+def atm_to_pa(p):
+    return p*101.325*1e3
+
+def psi_to_pa(p):
+    return p*6894.8
+
+# convert from Pa
+def pa_to_bar(p):
+    return p*1e-5
+
+def pa_to_mWS(p):
+    return p*1/(9.80665*1e3)
+
+def pa_to_torr(p):
+    return p/133.322
+
+def pa_to_atm(p):
+    return p*1/(101.325*1e3)
+
+ # converstion function
+
+def pa_to_psi(p):
+    return p/6894.8
+
+def pressure_conversion(pressure, input_unit = 'bar', target_unit = 'Pa'):
+    p = pressure
+    if input_unit.lower()   == 'bar':
+        p_pa = bar_to_pa(p)
+    elif input_unit.lower() == 'mws':
+        p_pa = mWS_to_pa(p)
+    elif input_unit.lower() == 'torr':
+        p_pa = torr_to_pa(p)
+    elif input_unit.lower() == 'atm':
+        p_pa = atm_to_pa(p)
+    elif input_unit.lower() == 'psi':
+        p_pa = psi_to_pa(p)
+    elif input_unit.lower() == 'pa':
+        p_pa = p
+    else:
+        raise Exception('Given input unit not recognised. \n Known units are: Pa, bar, mWs, Torr, atm, psi')
+    
+    if target_unit.lower()      == 'bar':
+        return pa_to_bar(p_pa), target_unit
+    elif target_unit.lower()    == 'mws':
+        return pa_to_mWS(p_pa), target_unit
+    elif target_unit.lower()    == 'torr':
+        return pa_to_torr(p_pa), target_unit
+    elif target_unit.lower()    == 'atm':
+        return pa_to_atm(p_pa), target_unit
+    elif target_unit.lower()    =='psi':
+        return pa_to_psi(p_pa), target_unit
+    elif target_unit.lower()    == 'pa':
+        return p_pa, target_unit
+    else:
+        raise Exception('Given target unit not recognised. \n Known units are: Pa, bar, mWs, Torr, atm, psi')
+
+# testing_pressure_conversion
+if __name__ == '__main__':
+    p = 1
+
+    unit_dict = ['Pa','Bar','Torr','Atm','MWS','psi']
+
+    for input_unit in unit_dict:
+        for target_unit in unit_dict:
+            converted_p = pressure_conversion(p,input_unit,target_unit)
+            print(input_unit,target_unit)
+            print(converted_p)

Druckrohrleitung/Druckstoß.ipynb

@@ -0,0 +1,166 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 37,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#imports\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 38,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#define constants\n",
+    "\n",
+    "g       = 9.81              # gravitational acceleration [m/s²]\n",
+    "\n",
+    "L       = 100              # length of pipeline [m]\n",
+    "rho     = 1000              # density of water [kg/m³]\n",
+    "D       = 1                 # pipe diameter \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      = 150              # number of time steps\n",
+    "f_coeff = 0.1              # lambda = 0.01 Friction loss coefficient [m]\n",
+    "c       = 400               # propagation velocity of the pressure wave\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 39,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 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_coeff*pl_vec/D*rho/2*v0**2)   # ref Wikipedia: Rohrreibungszahls\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 nn (at reservoir)\n",
+    "p_nn = np.zeros_like(t_vec)\n",
+    "v_nn = np.zeros_like(t_vec)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 40,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib qt\n",
+    "# time loop\n",
+    "\n",
+    "fig = plt.figure()\n",
+    "ax1 = fig.add_subplot(111)\n",
+    "lo1, = ax1.plot(pl_vec,np.full_like(pl_vec,p0),marker='.')\n",
+    "ax1.set_ylim([-20*p0,20*p0])\n",
+    "\n",
+    "for it in range(nt):\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_coeff*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_coeff*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_coeff*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_coeff*dt/(4*D)*(abs(v_old[i-1])*v_old[i-1]-abs(v_old[i+1])*v_old[i+1])\n",
+    "    \n",
+    "    lo1.set_xdata(pl_vec)\n",
+    "    lo1.set_ydata(p_new)\n",
+    "    ax1.set_title(str(t_vec[it]))\n",
+    "    fig.canvas.draw()\n",
+    "    plt.pause(0.05)\n",
+    "\n",
+    "    # store parameters of node 0 (at reservoir)\n",
+    "    p_nn[it] = p_old[-1]\n",
+    "    v_nn[it] = v_old[-1]\n",
+    "\n",
+    "    # prepare for next loop\n",
+    "    p_old = p_new\n",
+    "    v_old = v_new\n",
+    "\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 42,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[<matplotlib.lines.Line2D at 0x14c45f0bfd0>]"
+      ]
+     },
+     "execution_count": 42,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "plt.plot(v_nn)"
+   ]
+  }
+ ],
+ "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
+}