Georg_Brantegger/Python-DT_Slot_3
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'Dev'

Browse Source
This commit is contained in:
Brantegger Georg
2023-02-06 11:05:39 +01:00
parent 9f2c96088b 52c3a3c87d
commit bd2d17ee7c
156 changed files with 7126042 additions and 296 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
Druckrohrleitung/Druckrohrleitung_class_file.py
View File

@@ -1,13 +1,18 @@
# import modules for general use
import os # to import functions from other folders
import sys # to import functions from other folders
from logging import \
exception # to throw an exception when a specific condition is met
import numpy as np
#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 Druckrohrleitung_class:
# units
# make sure that units and display units are the same
@@ -37,6 +42,7 @@ class Druckrohrleitung_class:
g = 9.81 # m/s² gravitational acceleration
# init
# see docstring below
def __init__(self,total_length,diameter,pipeline_head,number_segments,Darcy_friction_factor,pw_vel,timestep,pressure_unit_disp,rho=1000):
"""
Creates a reservoir with given attributes in this order: \n
@@ -152,6 +158,7 @@ class Druckrohrleitung_class:
ss_v0 = np.full_like(self.x_vec,ss_flux/self.A)
# the static pressure is given by static state pressure of the reservoir, corrected for the hydraulic head of the pipe and friction losses
# dynamic pressure does not play a role, because it has the same influence on both sides of the equation (constant flow velocity) and therefore cancels out
ss_pressure = ss_pressure_res+(self.density*self.g*self.h_vec)-(self.f_D*self.x_vec/self.dia*self.density/2*ss_v0**2)
# set the initial conditions
@@ -160,6 +167,7 @@ class Druckrohrleitung_class:
# getter - return attributes
def get_info(self):
# prints out the info on the current state of the reservoir
new_line = '\n'
angle_deg = round(self.angle/np.pi*180,3)
@@ -180,8 +188,11 @@ class Druckrohrleitung_class:
f"Pressure wave vel. = {self.c:<10} {self.velocity_unit_disp} {new_line}"
f"Simulation timestep = {self.dt:<10} {self.time_unit_disp} {new_line}"
f"----------------------------- {new_line}"
f"Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object")
f"Velocity and pressure distribution are vectors and are accessible via the {new_line} \
get_current_velocity_distribution() and get_current_pressure_distribution() methods of the pipeline object. {new_line} \
See also get_lowest_XXX_per_node() and get_highest_XXX_per_node() methods.")
# print the info to console
print(print_str)
def get_current_pressure_distribution(self,disp_flag=False):
@@ -198,12 +209,14 @@ class Druckrohrleitung_class:
return self.v*self.A
def get_lowest_pressure_per_node(self,disp_flag=False):
# disp_flag if one wants to directly plot the return of this method
if disp_flag == True: # convert to pressure unit disp
return pressure_conversion(self.p_min,self.pressure_unit,self.pressure_unit_disp)
elif disp_flag == False: # stay in Pa
return self.p_min
def get_highest_pressure_per_node(self,disp_flag=False):
# disp_flag if one wants to directly plot the return of this method
if disp_flag == True: # convert to pressure unit disp
return pressure_conversion(self.p_max,self.pressure_unit,self.pressure_unit_disp)
elif disp_flag == False: # stay in Pa
@@ -229,8 +242,8 @@ class Druckrohrleitung_class:
return self.p0
def timestep_characteristic_method(self):
# use the method of characteristics to calculate the pressure and velocities at all nodes except the boundary ones
# they are set with the .set_boundary_conditions_next_timestep() method beforehand
# use the method of characteristics to calculate the pressure and velocities at all nodes except the boundary ones
# they are set with the .set_boundary_conditions_next_timestep() method beforehand
# constants for cleaner formula
nn = self.n_seg+1 # number of nodes
@@ -242,7 +255,7 @@ class Druckrohrleitung_class:
g = self.g # graviational acceleration
alpha = self.angle # pipeline angle
# Vectorize this loop?
# Vectorized loop see below
for i in range(1,nn-1):
self.v[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]) \
+dt*g*np.sin(alpha)-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])
@@ -263,8 +276,9 @@ class Druckrohrleitung_class:
self.v_old = self.v.copy()
def timestep_characteristic_method_vectorized(self):
# use the method of characteristics to calculate the pressure and velocities at all nodes except the boundary ones
# they are set with the .set_boundary_conditions_next_timestep() method beforehand
# faster then above
# use the method of characteristics to calculate the pressure and velocities at all nodes except the boundary ones
# they are set with the .set_boundary_conditions_next_timestep() method beforehand
# constants for cleaner formula
rho = self.density # density of liquid

142
Druckrohrleitung/Druckstoß Visualisierung.ipynb
View File

@@ -70,7 +70,7 @@
" # for general simulation\n",
"flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n",
"level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n",
"simTime_target = 10. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n",
"simTime_target = 3. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n",
"nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n",
"t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n"
]
@@ -79,6 +79,23 @@
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"61.1829727786757\n"
]
}
],
"source": [
"print(pressure_conversion(600000,'Pa','mWS'))"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# create objects\n",
@@ -94,18 +111,20 @@
},
{
"cell_type": "code",
"execution_count": 12,
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# initialization for timeloop\n",
"reservoir.set_influx = 0.\n",
"\n",
"level_vec = np.zeros_like(t_vec)\n",
"level_vec[0] = reservoir.get_current_level()\n",
"\n",
"# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n",
"v_old = pipe.get_current_velocity_distribution()\n",
"p_old = pipe.get_current_pressure_distribution()\n",
"v_old = pipe.get_current_velocity_distribution()\n",
"p_old = pipe.get_current_pressure_distribution()\n",
"p_0 = pipe.get_initial_pressure_distribution()\n",
"\n",
"# prepare the vectors in which the temporal evolution of the boundary conditions are stored\n",
" # keep in mind, that the velocity at the turbine and the pressure at the reservoir are set manually and\n",
@@ -135,30 +154,6 @@
"# v_boundary_tur[np.argmin(np.abs(t_vec-1)):] = 0"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x1efa21574f0>]"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"%matplotlib qt5\n",
"fig = plt.figure()\n",
"plt.plot(v_trans)\n",
"fig = plt.figure()\n",
"plt.plot(t_vec,v_boundary_tur)"
]
},
{
"cell_type": "code",
"execution_count": 14,
@@ -166,21 +161,35 @@
"outputs": [],
"source": [
"%matplotlib qt5\n",
"fig1,axs1 = plt.subplots(2,1)\n",
"# Time loop\n",
"\n",
"# create a figure and subplots to display the velocity and pressure distribution across the pipeline in each pipeline step\n",
"fig1,axs1 = plt.subplots(3,1, figsize=(16,9))\n",
"fig1.suptitle(str(0) +' s / '+str(round(t_vec[-1],2)) + ' s' )\n",
"axs1[0].set_title('Pressure distribution in pipeline')\n",
"axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs1[0].set_ylabel(r'$p$ [mWS]')\n",
"axs1[0].autoscale()\n",
"lo_00, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,'Pa',pUnit_conv),marker='.')\n",
"\n",
"axs1[1].set_title('Velocity distribution in pipeline')\n",
"axs1[0].set_ylabel(r'$p$ ['+pUnit_conv+']')\n",
"axs1[0].set_ylim([-2,200])\n",
"axs1[1].set_title('Pressure distribution in pipeline \\n Difference to t=0')\n",
"axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs1[1].set_ylabel(r'$v$ [m/s]')\n",
"lo_01, = axs1[1].plot(Pip_x_vec,v_old,marker='.')\n",
"# axs1[1].autoscale()\n",
"axs1[1].set_ylim([-1.5,1.5])\n",
"axs1[1].set_ylabel(r'$p$ ['+pUnit_conv+']')\n",
"axs1[1].set_ylim([-76,76])\n",
"axs1[2].set_title('Flux distribution in pipeline')\n",
"axs1[2].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"axs1[2].set_ylabel(r'$Q$ [$\\mathrm{m}^3 / \\mathrm{s}$]')\n",
"axs1[2].set_ylim([-1.5,1.5])\n",
"lo_0, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,pUnit_calc, pUnit_conv),marker='.')\n",
"lo_0min, = axs1[0].plot(Pip_x_vec,pressure_conversion(pipe.get_lowest_pressure_per_node(),pUnit_calc,pUnit_conv),c='red')\n",
"lo_0max, = axs1[0].plot(Pip_x_vec,pressure_conversion(pipe.get_highest_pressure_per_node(),pUnit_calc,pUnit_conv),c='red')\n",
"lo_1, = axs1[1].plot(Pip_x_vec,pressure_conversion(p_old-p_0,pUnit_calc, pUnit_conv),marker='.')\n",
"lo_1min, = axs1[1].plot(Pip_x_vec,pressure_conversion(pipe.get_lowest_pressure_per_node()-p_0,pUnit_calc,pUnit_conv),c='red')\n",
"lo_1max, = axs1[1].plot(Pip_x_vec,pressure_conversion(pipe.get_highest_pressure_per_node()-p_0,pUnit_calc,pUnit_conv),c='red')\n",
"lo_2, = axs1[1].plot(Pip_x_vec,v_old,marker='.')\n",
"lo_2min, = axs1[2].plot(Pip_x_vec,pipe.get_lowest_velocity_per_node(),c='red')\n",
"lo_2max, = axs1[2].plot(Pip_x_vec,pipe.get_highest_velocity_per_node(),c='red')\n",
"\n",
"fig1.tight_layout()\n",
"fig1.show()\n",
"plt.pause(1)"
]
},
@@ -219,17 +228,38 @@
" # plot some stuff\n",
" if it_pipe%100 == 0:\n",
" # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n",
" lo_00.remove()\n",
" lo_01.remove()\n",
" # lo_02.remove()\n",
" lo_0.remove()\n",
" lo_0min.remove()\n",
" lo_0max.remove()\n",
" lo_1.remove()\n",
" lo_1min.remove()\n",
" lo_1max.remove()\n",
" lo_2.remove()\n",
" lo_2min.remove()\n",
" lo_2max.remove()\n",
" # plot new pressure and velocity distribution in the pipeline\n",
" lo_00, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,'Pa', pUnit_conv),marker='.',c='blue')\n",
" lo_01, = axs1[1].plot(Pip_x_vec,v_old,marker='.',c='blue')\n",
" \n",
" fig1.suptitle(str(round(t_vec[it_pipe],2)) + '/' + str(round(t_vec[-1],2)))\n",
" lo_0, = axs1[0].plot(Pip_x_vec,pressure_conversion(pipe.get_current_pressure_distribution(),pUnit_calc,pUnit_conv),marker='.',c='blue')\n",
" lo_0min, = axs1[0].plot(Pip_x_vec,pressure_conversion(pipe.get_lowest_pressure_per_node(),pUnit_calc,pUnit_conv),c='red')\n",
" lo_0max, = axs1[0].plot(Pip_x_vec,pressure_conversion(pipe.get_highest_pressure_per_node(),pUnit_calc,pUnit_conv),c='red') \n",
" lo_1, = axs1[1].plot(Pip_x_vec,pressure_conversion(pipe.get_current_pressure_distribution()-p_0,pUnit_calc,pUnit_conv),marker='.',c='blue')\n",
" lo_1min, = axs1[1].plot(Pip_x_vec,pressure_conversion(pipe.get_lowest_pressure_per_node()-p_0,pUnit_calc,pUnit_conv),c='red')\n",
" lo_1max, = axs1[1].plot(Pip_x_vec,pressure_conversion(pipe.get_highest_pressure_per_node()-p_0,pUnit_calc,pUnit_conv),c='red')\n",
" lo_2, = axs1[2].plot(Pip_x_vec,pipe.get_current_flux_distribution(),marker='.',c='blue')\n",
" lo_2min, = axs1[2].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n",
" lo_2max, = axs1[2].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n",
" fig1.suptitle(str(round(t_vec[it_pipe]-1,2))+ ' s / '+str(round(t_vec[-1]-1,2)) + ' s' )\n",
" fig1.canvas.draw()\n",
" fig1.tight_layout()\n",
" plt.pause(0.000001)"
" fig1.show()\n",
" # if int(it_pipe/100) < 10:\n",
" # figname = 'GIF Plots\\ GIF00'+str(int(it_pipe/100))+'.png'\n",
" # elif int(it_pipe/100) < 100:\n",
" # figname = 'GIF Plots\\ GIF0'+str(int(it_pipe/100))+'.png'\n",
" # else:\n",
" # figname = 'GIF Plots\\ GIF'+str(int(it_pipe/100))+'.png'\n",
" # print(figname)\n",
" # fig1.savefig(fname=figname)\n",
" plt.pause(0.000001) "
]
},
{
@@ -245,11 +275,11 @@
"axs2[0,0].set_ylabel(r'$p$ [mWS]')\n",
"axs2[0,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_res,pUnit_calc,pUnit_conv)),1.1*np.max(pressure_conversion(p_boundary_res,pUnit_calc,pUnit_conv))])\n",
"\n",
"axs2[1,1].set_title('Velocity Reservoir')\n",
"axs2[1,1].plot(t_vec,v_boundary_res)\n",
"axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
"axs2[1,1].set_ylim([-1.1*np.max(v_boundary_res),1.1*np.max(v_boundary_res)])\n",
"axs2[1,0].set_title('Velocity Reservoir')\n",
"axs2[1,0].plot(t_vec,v_boundary_res)\n",
"axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[1,0].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
"axs2[1,0].set_ylim([-1.1*np.max(v_boundary_res),1.1*np.max(v_boundary_res)])\n",
"\n",
"axs2[0,1].set_title('Pressure Turbine')\n",
"axs2[0,1].plot(t_vec,pressure_conversion(p_boundary_tur,pUnit_calc,pUnit_conv))\n",
@@ -257,11 +287,11 @@
"axs2[0,1].set_ylabel(r'$p$ [mWS]')\n",
"axs2[0,1].set_ylim([0.9*np.min(pressure_conversion(p_boundary_tur,pUnit_calc,pUnit_conv)),1.1*np.max(pressure_conversion(p_boundary_tur,pUnit_calc,pUnit_conv))])\n",
"\n",
"axs2[1,0].set_title('Velocity Turbine')\n",
"axs2[1,0].plot(t_vec,v_boundary_tur)\n",
"axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[1,0].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
"axs2[1,0].set_ylim([-0.1,1.05*np.max(v_boundary_tur)])\n",
"axs2[1,1].set_title('Velocity Turbine')\n",
"axs2[1,1].plot(t_vec,v_boundary_tur)\n",
"axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
"axs2[1,1].set_ylim([-0.1,1.05*np.max(v_boundary_tur)])\n",
"\n",
"fig2.tight_layout()\n",
"plt.show()"