added first try for turbine flux based on pressure and

LA opening

Turbinen/Turbinen_class_file.py

@@ -1,7 +1,4 @@
-from matplotlib.pyplot import fill
 import numpy as np
-from scipy.interpolate import interp2d
-
 #importing pressure conversion function
 import sys
 import os
@@ -10,21 +7,28 @@ parent = os.path.dirname(current)
 sys.path.append(parent)
 from functions.pressure_conversion import pressure_conversion
 
-class Francis_turbine_class:
+class Francis_Turbine:
-    def __init__(self,CSV_name='Durchflusskennlinie.csv'):
+    def __init__(self, Q_nenn,p_nenn):
-        self.raw_csv = np.genfromtxt(CSV_name,delimiter=',')
+        self.Q_n    = Q_nenn
-        
+        self.p_n    = p_nenn
-    def extract_csv(self,CSV_pressure_unit='bar'):
+        self.LA_n   = 1. # 100%
-        self.raw_ps_vec,_    = pressure_conversion(self.raw_csv[0,1:],CSV_pressure_unit,'Pa')
+        h,_ = pressure_conversion(p_nenn,'Pa','MWs')
-        self.raw_LA_vec     = self.raw_csv[1:,0]
+        self.A      = Q_nenn/(np.sqrt(2*9.81*h)*0.98)
-        self.raw_Qs_mat      = self.raw_csv[1:,1:]
-
-    def get_Q_fun(self):
-        Q_fun = interp2d(self.raw_ps_vec,self.raw_LA_vec,self.raw_Qs_mat,bounds_error=False,fill_value=None)
-        return Q_fun
-
-
 
+    def set_LA(self,LA):
+        self.LA = LA
 
+    def get_Q(self,p):
+        self.Q = self.Q_n*(self.LA/self.LA_n)*np.sqrt(p/self.p_n)
+        return self.Q
 
+    def set_closing_time(self,t_closing):
+        self.t_c = t_closing
+        self.d_LA_max_dt = 1/t_closing
 
+    def change_LA(self,LA_soll,timestep):
+        LA_diff = self.LA-LA_soll
+        LA_diff_max = self.d_LA_max_dt*timestep
+        if abs(LA_diff) > LA_diff_max:
+            LA_diff = np.sign(LA_diff)*LA_diff_max
+        self.LA = self.LA-LA_diff

Turbinen/old/Turbinen_class_file.py

@@ -0,0 +1,33 @@
+from matplotlib.pyplot import fill
+import numpy as np
+from scipy.interpolate import interp2d
+
+#importing pressure conversion function
+import sys
+import os
+current = os.path.dirname(os.path.realpath(__file__))
+parent = os.path.dirname(current)
+sys.path.append(parent)
+from functions.pressure_conversion import pressure_conversion
+
+class Francis_turbine_class:
+    def __init__(self,CSV_name='Durchflusskennlinie.csv'):
+        csv = np.genfromtxt(CSV_name,delimiter=',')
+        n_rows,_ = np.shape(csv)
+        self.raw_csv = np.append(csv,np.zeros([n_rows,1]),axis = 1)
+        
+    def extract_csv(self,CSV_pressure_unit='bar'):
+        ps_vec,_            = pressure_conversion(self.raw_csv[0,1:],CSV_pressure_unit,'Pa')
+        self.raw_ps_vec     = np.flip(ps_vec)
+        self.raw_LA_vec     = self.raw_csv[1:,0]
+        self.raw_Qs_mat     = np.fliplr(self.raw_csv[1:,1:])/1000. # convert from l/s to m³/s
+
+    def get_Q_fun(self):
+        Q_fun = interp2d(self.raw_ps_vec,self.raw_LA_vec,self.raw_Qs_mat,bounds_error=False,fill_value=None)
+        return Q_fun
+
+
+
+
+
+

Untertweng.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 56,
+   "execution_count": 46,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -11,41 +11,46 @@
     "\n",
     "from functions.pressure_conversion import pressure_conversion\n",
     "from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class\n",
-    "from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class"
+    "from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class\n",
+    "from Turbinen.Turbinen_class_file import Francis_Turbine"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 57,
+   "execution_count": 47,
    "metadata": {},
    "outputs": [],
    "source": [
     "#define constants\n",
     "\n",
+    "#Turbine\n",
+    "Q_nenn = 0.85\n",
+    "p_nenn,_ = pressure_conversion(10.6,'bar','Pa')\n",
+    "\n",
     "# physics\n",
     "g               = 9.81                          # gravitational acceleration [m/s²]\n",
     "rho             = 1000.                         # density of water [kg/m³]\n",
     "\n",
     "# pipeline\n",
-    "L               = 1000.                         # length of pipeline [m]\n",
+    "L               = 535.+478.                        # length of pipeline [m]\n",
-    "D               = 1.                            # pipe diameter [m]\n",
+    "D               = 0.9                            # pipe diameter [m]\n",
     "A_pipe          = D**2/4*np.pi                  # pipeline area\n",
-    "h_pipe          = 200                           # hydraulic head without reservoir [m] \n",
+    "h_pipe          = 105                           # hydraulic head without reservoir [m] \n",
     "alpha           = np.arcsin(h_pipe/L)           # Höhenwinkel der Druckrohrleitung \n",
     "n               = 50                           # number of pipe segments in discretization\n",
     "# consider replacing Q0 with a vector be be more flexible in initial conditions\n",
-    "Q0              = 2.                            # initial flow in whole pipe [m³/s]\n",
+    "Q0              = Q_nenn                            # initial flow in whole pipe [m³/s]\n",
     "v0              = Q0/A_pipe                     # initial flow velocity [m/s]\n",
-    "f_D             = 0.01                           # Darcy friction factor\n",
+    "f_D             = 0.014                           # Darcy friction factor\n",
-    "c               = 400.                          # propagation velocity of the pressure wave [m/s]\n",
+    "c               = 500.                          # propagation velocity of the pressure wave [m/s]\n",
     "# consider prescribing a total simulation time and deducting the number of timesteps from that\n",
-    "nt              = 500                           # number of time steps after initial conditions\n",
+    "nt              = 2000                           # number of time steps after initial conditions\n",
     "\n",
     "# derivatives of the pipeline constants\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",
-    "initial_level   = 20.                           # water level in upstream reservoir [m]\n",
+    "initial_level   = 8.                           # water level in upstream reservoir [m]\n",
     "p0              = rho*g*initial_level-v0**2*rho/2\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+1)*dt          # time vector\n",
@@ -61,7 +66,7 @@
     "initial_pipeline_pressure   = p0        # Initial condition for the static pipeline pressure at the reservoir (= hydrostatic pressure - dynamic pressure) \n",
     "initial_pressure_unit       = 'Pa'      # DO NOT CHANGE! for pressure conversion in print statements and plot labels \n",
     "conversion_pressure_unit    = 'bar'      # for pressure conversion in print statements and plot labels\n",
-    "area_base                   = 20.       # total base are of the cuboid reservoir [m²]   \n",
+    "area_base                   = 74.       # total base are of the cuboid reservoir [m²]   \n",
     "area_outflux                = A_pipe    # outlfux area of the reservoir, given by pipeline area [m²]\n",
     "critical_level_low          = 0.        # for yet-to-be-implemented warnings[m]\n",
     "critical_level_high         = np.inf    # for yet-to-be-implemented warnings[m]\n",
@@ -69,6 +74,7 @@
     "# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n",
     "nt_eRK4                     = 1000      # number of simulation steps of reservoir in between timesteps of pipeline              \n",
     "simulation_timestep         = dt/nt_eRK4\n",
+    "\n",
     "\n"
    ]
   },
@@ -91,7 +97,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 58,
+   "execution_count": 48,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -109,6 +115,11 @@
     "pipe.set_initial_pressure(p_init,initial_pressure_unit,conversion_pressure_unit)\n",
     "pipe.set_initial_flow_velocity(v_init)\n",
     "\n",
+    "\n",
+    "T1 = Francis_Turbine(Q_nenn,p_nenn)\n",
+    "T1.set_LA(1.)\n",
+    "T1.set_closing_time(30)\n",
+    "\n",
     "# display the attributes of the created reservoir and pipeline object\n",
     "# V.get_info(full=True)\n",
     "# pipe.get_info()"
@@ -116,7 +127,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 59,
+   "execution_count": 49,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -139,21 +150,21 @@
     "level_vec       = np.full(nt+1,initial_level)  # level at the end of each pipeline timestep\n",
     "level_vec_2     = np.empty([nt_eRK4])   # level throughout each reservoir timestep-used for plotting and overwritten afterwards\n",
     "\n",
-    "# set the boudary conditions for the first timestep\n",
+    "# set the boundary conditions for the first timestep\n",
     "v_boundary_res[0]   = v_old[0]\n",
     "v_boundary_tur[0]   = v_old[-1] \n",
-    "v_boundary_tur[1:]  = 0                                     # instantaneous closing\n",
+    "p_boundary_res[0]   = p_old[0]\n",
-    "# v_boundary_tur[0:20]    = np.linspace(v_old[-1],0,20)    # overwrite for finite closing time - linear case\n",
+    "p_boundary_tur[0]   = p_old[-1]\n",
-    "# const = int(np.min([100,round(nt/1.1)]))\n",
+    "\n",
-    "# v_boundary_tur[0:const] = v_old[1]*np.cos(t_vec[0:const]*2*np.pi/5)**2\n",
+    "LA_soll_vec = np.zeros_like(t_vec)\n",
-    "p_boundary_res[0]       = p_old[0]\n",
+    "LA_soll_vec[0] = 1\n",
-    "p_boundary_tur[0]       = p_old[-1]\n",
+    "LA_soll_vec[1000:] = 1\n",
     "\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 60,
+   "execution_count": 50,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -186,7 +197,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 61,
+   "execution_count": 51,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -213,6 +224,9 @@
     "    v_boundary_res[it_pipe] = v_old[1]+1/(rho*c)*(p_boundary_res[it_pipe]-p_old[1])-f_D*dt/(2*D)*abs(v_old[1])*v_old[1] \\\n",
     "        +dt*g*np.sin(alpha)\n",
     "\n",
+    "    T1.change_LA(LA_soll_vec[it_pipe],dt)\n",
+    "    v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_Q(p_old[-1])\n",
+    "\n",
     "    # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n",
     "    pipe.set_boundary_conditions_next_timestep(v_boundary_res[it_pipe],p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n",
     "    p_boundary_tur[it_pipe] = pipe.p_boundary_tur\n",
@@ -245,7 +259,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 62,
+   "execution_count": 52,
    "metadata": {},
    "outputs": [],
    "source": [

untertweng.txt

@@ -18,3 +18,5 @@ c = 500 m/s
 
 Q_0 = 100%*0.75+30%*0.75
 Q_extrem = 30%*0.75
+
+Q = LA*Q_nenn*sqrt(H/H_n)