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

code cleanup and commentary

Browse Source
This commit is contained in:
Brantegger Georg
2022-07-06 10:29:22 +02:00
parent b03bb43c63
commit c81c0ab142
6 changed files with 188 additions and 156 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
Ausgleichsbecken/Ausgleichsbecken_class_file.py
View File

@@ -66,13 +66,14 @@ class Ausgleichsbecken_class:
# getter # getter
def get_info(self, full = False): def get_info(self, full = False):
new_line = '\n' new_line = '\n'
if full == True: if full == True:
# :<10 pads the self.value to be 10 characters wide # :<10 pads the self.value to be 10 characters wide
print_str = (f"The cuboid reservoir has the following attributes: {new_line}" print_str = (f"The cuboid reservoir has the following attributes: {new_line}"
f"----------------------------- {new_line}" f"----------------------------- {new_line}"
f"Base area = {self.area:<10} {self.area_unit_print} {new_line}" f"Base area = {self.area:<10} {self.area_unit_print} {new_line}"
f"Outflux area = {self.area_outflux:<10} {self.area_outflux_unit_print} {new_line}" f"Outflux area = {round(self.area_outflux,3):<10} {self.area_outflux_unit_print} {new_line}"
f"Current level = {self.level:<10} {self.level_unit_print}{new_line}" f"Current level = {self.level:<10} {self.level_unit_print}{new_line}"
f"Critical level low = {self.level_min:<10} {self.level_unit_print} {new_line}" f"Critical level low = {self.level_min:<10} {self.level_unit_print} {new_line}"
f"Critical level high = {self.level_max:<10} {self.level_unit_print} {new_line}" f"Critical level high = {self.level_max:<10} {self.level_unit_print} {new_line}"

30
Ausgleichsbecken/Ausgleichsbecken_functions.py
View File

@@ -1,30 +0,0 @@
import numpy as np
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, p0)+2*FODE_function(Y2, h_hs, alpha, p_hs)+ \
2*FODE_function(Y3, h_hs, alpha, p_hs)+ FODE_function(Y4, h, alpha, p0))

23
Druckrohrleitung/Druckrohrleitung_class_file.py
View File

@@ -9,6 +9,8 @@ sys.path.append(parent)
from functions.pressure_conversion import pressure_conversion from functions.pressure_conversion import pressure_conversion
class Druckrohrleitung_class: class Druckrohrleitung_class:
# units # units
acceleration_unit = r'$\mathrm{m}/\mathrm{s}^2$' acceleration_unit = r'$\mathrm{m}/\mathrm{s}^2$'
@@ -109,7 +111,9 @@ class Druckrohrleitung_class:
# getter # getter
def get_info(self): def get_info(self):
new_line = '\n' new_line = '\n'
angle_deg = round(self.angle/np.pi*180,3)
# :<10 pads the self.value to be 10 characters wide # :<10 pads the self.value to be 10 characters wide
print_str = (f"The pipeline has the following attributes: {new_line}" print_str = (f"The pipeline has the following attributes: {new_line}"
@@ -118,13 +122,15 @@ class Druckrohrleitung_class:
f"Diameter = {self.dia:<10} {self.length_unit_print} {new_line}" f"Diameter = {self.dia:<10} {self.length_unit_print} {new_line}"
f"Number of segments = {self.n_seg:<10} {new_line}" f"Number of segments = {self.n_seg:<10} {new_line}"
f"Number of nodes = {self.n_seg+1:<10} {new_line}" f"Number of nodes = {self.n_seg+1:<10} {new_line}"
f"Length per segments = {self.dx:<10} {self.length_unit_print} {new_line}" f"Length per segments = {self.dx:<10} {self.length_unit_print} {new_line}"
f"Pipeline angle = {round(self.angle,3):<10} {self.angle_unit_print} {new_line}" f"Pipeline angle = {round(self.angle,3):<10} {self.angle_unit_print} {new_line}"
f"Pipeline angle = {angle_deg}° {new_line}"
f"Darcy friction factor = {self.f_D:<10} {new_line}" f"Darcy friction factor = {self.f_D:<10} {new_line}"
f"Density of liquid = {self.density:<10} {self.density_unit_print} {new_line}" f"Density of liquid = {self.density:<10} {self.density_unit_print} {new_line}"
f"Pressure wave vel. = {self.c:<10} {self.velocity_unit_print} {new_line}" f"Pressure wave vel. = {self.c:<10} {self.velocity_unit_print} {new_line}"
f"Simulation timestep = {self.dt:<10} {self.time_unit_print } {new_line}" f"Simulation timestep = {self.dt:<10} {self.time_unit_print} {new_line}"
f"Number of timesteps = {self.nt:<10} {new_line}" f"Number of timesteps = {self.nt:<10} {new_line}"
f"Total simulation time = {self.nt*self.dt:<10} {self.time_unit_print} {new_line}"
f"----------------------------- {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 by the .v and .p attribute of the pipeline object")
@@ -162,14 +168,3 @@ class Druckrohrleitung_class:
self.p_old = self.p.copy() self.p_old = self.p.copy()
self.v_old = self.v.copy() self.v_old = self.v.copy()

0
Druckrohrleitung/Druckrohrleitung_ETH_class_file.py → Druckrohrleitung/old/Druckrohrleitung_ETH_class_file.py
View File

0
Druckrohrleitung/Druckstoß_ETH.ipynb → Druckrohrleitung/old/Druckstoß_ETH.ipynb
View File

288
combine_pipeline_and_reservoir.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 5, "execution_count": 6,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -16,7 +16,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 6, "execution_count": 7,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -24,47 +24,48 @@
"\n", "\n",
"# physics\n", "# physics\n",
"g = 9.81 # gravitational acceleration [m/s²]\n", "g = 9.81 # gravitational acceleration [m/s²]\n",
"rho = 1000. # density of water [kg/m³]\n", "rho = 1000. # density of water [kg/m³]\n",
"\n", "\n",
"# pipeline\n", "# pipeline\n",
"L = 1000. # length of pipeline [m]\n", "L = 1000. # length of pipeline [m]\n",
"D = 1. # pipe diameter [m]\n", "D = 1. # pipe diameter [m]\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",
"A_pipe = D**2/4*np.pi # pipeline area\n", "A_pipe = D**2/4*np.pi # pipeline area\n",
"v0 = Q0/A_pipe # initial flow velocity [m/s]\n", "h_pipe = 200 # hydraulic head without reservoir [m] \n",
"h_res = 20. # water level in upstream reservoir [m]\n",
"n = 10 # number of pipe segments in discretization\n",
"nt = 10000 # 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",
"h_pipe = 300 # hydraulic head without reservoir [m] \n",
"alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n", "alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n",
"n = 10 # 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",
"v0 = Q0/A_pipe # initial flow velocity [m/s]\n",
"f_D = 0.1 # Darcy friction factor\n",
"c = 400. # 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",
"\n", "\n",
"# derivatives of the pipeline constants\n", "# derivatives of the pipeline constants\n",
"p0 = rho*g*h_res-v0**2*rho/2\n", "dx = L/n # length of each pipe segment\n",
"dx = L/n # length of each pipe segment\n", "dt = dx/c # timestep according to method of characterisitics\n",
"dt = dx/c # timestep according to method of characterisitics\n", "nn = n+1 # number of nodes\n",
"nn = n+1 # number of nodes\n", "h_res = 20. # water level in upstream reservoir [m]\n",
"pl_vec = np.arange(0,nn*dx,dx) # pl = pipe-length. position of the nodes on the pipeline\n", "p0 = rho*g*h_res-v0**2*rho/2\n",
"t_vec = np.arange(0,nt*dt,dt) # time vector\n", "pl_vec = np.arange(0,nn*dx,dx) # pl = pipe-length. position of the nodes on the pipeline\n",
"h_vec = np.arange(0,h_pipe+h_pipe/n,h_pipe/n) # hydraulic head of pipeline at each node\n", "t_vec = np.arange(0,nt*dt,dt) # time vector\n",
"\n", "h_vec = np.arange(0,n+1)*h_pipe/n # hydraulic head of pipeline at each node np.arange(0,0) does not yield the intended result\n",
"v_init = np.full(nn,Q0/(D**2/4*np.pi))\n", "v_init = np.full(nn,Q0/(D**2/4*np.pi)) # initial velocity distribution in pipeline\n",
"p_init = (rho*g*(h_res+h_vec)-v_init**2*rho/2)-(f_D*pl_vec/D*rho/2*v_init**2) # ref Wikipedia: Darcy Weisbach\n", "p_init = (rho*g*(h_res+h_vec)-v_init**2*rho/2)-(f_D*pl_vec/D*rho/2*v_init**2) # ref Wikipedia: Darcy Weisbach\n",
"\n", "\n",
"\n", "\n",
"# reservoir\n", "# reservoir\n",
"initial_level = h_res # m\n", "initial_level = h_res # water level in upstream reservoir [m]\n",
"initial_influx = 0. # m³/s\n", "# replace influx by vector\n",
"initial_outflux = Q0 # m³/s\n", "initial_influx = 0. # initial influx of volume to the reservoir [m³/s]\n",
"initial_pipeline_pressure = p0 # Pa \n", "initial_outflux = Q0 # initial outflux of volume from the reservoir to the pipeline [m³/s]\n",
"initial_pressure_unit = 'Pa'\n", "initial_pipeline_pressure = p0 # Initial condition for the static pipeline pressure at the reservoir (= hydrostatic pressure - dynamic pressure) \n",
"conversion_pressure_unit = 'Pa'\n", "initial_pressure_unit = 'Pa' # for pressure conversion in print statements and plot labels\n",
"area_base = 5. # m² really large base are to ensure level never becomes < 0\n", "conversion_pressure_unit = 'Pa' # for pressure conversion in print statements and plot labels\n",
"area_outflux = A_pipe # m²\n", "area_base = 20. # total base are of the cuboid reservoir [m²] \n",
"critical_level_low = 0. # m\n", "area_outflux = A_pipe # outlfux area of the reservoir, given by pipeline area [m²]\n",
"critical_level_high = np.inf # 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",
"\n", "\n",
"# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n", "# 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", "nt_eRK4 = 1000 # number of simulation steps of reservoir in between timesteps of pipeline \n",
@@ -73,27 +74,64 @@
] ]
}, },
{ {
"cell_type": "code", "cell_type": "markdown",
"execution_count": 7,
"metadata": {}, "metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(3.6368236494728476, 'mWS')\n"
]
}
],
"source": [ "source": [
"print(pressure_conversion(-np.sum((-v_init**2*rho/2)),'Pa','mWS'))" "#### Ideas for checks after constant definitions: \n",
"\n",
"- Check that the initial pressure is not negative:\n",
" - may happen, if there is too little hydraulic head to create the initial flow conditions with the given friction\n",
"<br>\n",
"<br>\n",
"- stupidity checks?\n",
" - area > area_outflux ?\n",
" - propable ranges for parameters?\n",
" - angle and height/length fit together?\n",
" "
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 8, "execution_count": 8,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The cuboid reservoir has the following attributes: \n",
"----------------------------- \n",
"Base area = 20.0 m² \n",
"Outflux area = 0.785 m² \n",
"Current level = 20.0 m\n",
"Critical level low = 0.0 m \n",
"Critical level high = inf m \n",
"Volume in reservoir = 400.0 m³ \n",
"Current influx = 0.0 m³/s \n",
"Current outflux = 2.0 m³/s \n",
"Simulation timestep = 0.00025 s \n",
"----------------------------- \n",
"\n",
"The pipeline has the following attributes: \n",
"----------------------------- \n",
"Length = 1000.0 m \n",
"Diameter = 1.0 m \n",
"Number of segments = 10 \n",
"Number of nodes = 11 \n",
"Length per segments = 100.0 m \n",
"Pipeline angle = 0.201 rad \n",
"Pipeline angle = 11.537° \n",
"Darcy friction factor = 0.1 \n",
"Density of liquid = 1000 kg/m³ \n",
"Pressure wave vel. = 400.0 m/s \n",
"Simulation timestep = 0.25 s \n",
"Number of timesteps = 500 \n",
"Total simulation time = 125.0 s \n",
"----------------------------- \n",
"Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object\n"
]
}
],
"source": [ "source": [
"# create objects\n", "# create objects\n",
"\n", "\n",
@@ -103,11 +141,16 @@
"V.set_outflux(initial_outflux)\n", "V.set_outflux(initial_outflux)\n",
"V.pressure, V.pressure_unit = pressure_conversion(initial_pipeline_pressure,input_unit = initial_pressure_unit, target_unit = conversion_pressure_unit)\n", "V.pressure, V.pressure_unit = pressure_conversion(initial_pipeline_pressure,input_unit = initial_pressure_unit, target_unit = conversion_pressure_unit)\n",
"\n", "\n",
"\n",
"pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n", "pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n",
"pipe.set_pressure_propagation_velocity(c)\n", "pipe.set_pressure_propagation_velocity(c)\n",
"pipe.set_number_of_timesteps(nt)\n", "pipe.set_number_of_timesteps(nt)\n",
"pipe.set_initial_pressure(p_init)\n", "pipe.set_initial_pressure(p_init)\n",
"pipe.set_initial_flow_velocity(v_init)" "pipe.set_initial_flow_velocity(v_init)\n",
"\n",
"# display the attributes of the created reservoir and pipeline object\n",
"V.get_info(full=True)\n",
"pipe.get_info()"
] ]
}, },
{ {
@@ -118,26 +161,33 @@
"source": [ "source": [
"# initialization for timeloop\n", "# initialization for timeloop\n",
"\n", "\n",
"# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n",
"v_old = v_init.copy()\n", "v_old = v_init.copy()\n",
"p_old = p_init.copy()\n", "p_old = p_init.copy()\n",
"\n", "\n",
"#vectors to store boundary conditions\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",
" # through the time evolution of the reservoir respectively \n",
" # the pressure at the turbine and the velocity at the reservoir are calculated from the method of characteristics\n",
"v_boundary_res = np.empty_like(t_vec)\n", "v_boundary_res = np.empty_like(t_vec)\n",
"v_boundary_tur = np.empty_like(t_vec)\n", "v_boundary_tur = np.empty_like(t_vec)\n",
"p_boundary_res = np.empty_like(t_vec)\n", "p_boundary_res = np.empty_like(t_vec)\n",
"p_boundary_tur = np.empty_like(t_vec)\n", "p_boundary_tur = np.empty_like(t_vec)\n",
"level_vec = np.empty_like(t_vec)\n",
"level_vec_2 = np.full([nt_eRK4],initial_level)\n",
"\n", "\n",
"# prepare the vectors that store the temporal evolution of the level in the reservoir\n",
"level_vec = np.full_like(t_vec,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",
"v_boundary_res[0] = v_old[0]\n", "v_boundary_res[0] = v_old[0]\n",
"v_boundary_tur[0] = v_old[-1] # instantaneous closing\n", "v_boundary_tur[0] = v_old[-1] \n",
"# v_boundary_tur[1:] = 0\n", "v_boundary_tur[1:] = 0 # instantaneous closing\n",
"v_boundary_tur[0:1000] = np.linspace(v_old[-1],0,1000) # finite closing time - linear case\n", "# v_boundary_tur[0:20] = np.linspace(v_old[-1],0,20) # overwrite for finite closing time - linear case\n",
"const = int(np.min([100,round(nt/1.1)]))\n",
"v_boundary_tur[0:const] = v_old[1]*np.cos(t_vec[0:const]*2*np.pi/5)**2\n",
"p_boundary_res[0] = p_old[0]\n", "p_boundary_res[0] = p_old[0]\n",
"p_boundary_tur[0] = p_old[-1]\n", "p_boundary_tur[0] = p_old[-1]\n",
"level_vec[0] = initial_level\n", "\n"
"\n",
"v_boundary_tur[1:] = 0 # instantaneous closing"
] ]
}, },
{ {
@@ -149,62 +199,73 @@
"%matplotlib qt5\n", "%matplotlib qt5\n",
"# time loop\n", "# time loop\n",
"\n", "\n",
"\n", "# create a figure and subplots to display the velocity and pressure distribution across the pipeline in each pipeline step\n",
"# fig2,axs2 = plt.subplots(3,1)\n", "fig1,axs1 = plt.subplots(2,1)\n",
"# axs2[0].set_title('Pressure distribution in pipeline')\n", "axs1[0].set_title('Pressure distribution in pipeline')\n",
"# axs2[1].set_title('Velocity distribution in pipeline')\n", "axs1[1].set_title('Velocity distribution in pipeline')\n",
"# axs2[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"# axs2[0].set_ylabel(r'$p$ [mWS]')\n", "axs1[0].set_ylabel(r'$p$ [mWS]')\n",
"# axs2[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
"# axs2[1].set_ylabel(r'$p$ [mWS]')\n", "axs1[1].set_ylabel(r'$v$ [$\\mathrm{m} / \\mathrm{s}$]')\n",
"# lo_00, = axs2[0].plot(pl_vec,pressure_conversion(pipe.p_old,'Pa','mWS')[0],marker='.')\n", "lo_00, = axs1[0].plot(pl_vec,pressure_conversion(pipe.p_old,'Pa','mWS')[0],marker='.')\n",
"# lo_01, = axs2[1].plot(pl_vec,pipe.v_old,marker='.')\n", "lo_01, = axs1[1].plot(pl_vec,pipe.v_old,marker='.')\n",
"# lo_02, = axs2[2].plot(level_vec_2)\n", "axs1[0].autoscale()\n",
"# axs2[0].autoscale()\n", "axs1[1].autoscale()\n",
"# axs2[1].autoscale()\n", "# displaying the reservoir level within each pipeline timestep\n",
"# axs2[2].autoscale()\n", "# axs1[2].set_title('Level reservoir')\n",
"# fig2.tight_layout()\n", "# axs1[2].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"# axs1[2].set_ylabel(r'$h$ [m]')\n",
"# lo_02, = axs1[2].plot(level_vec_2)\n",
"# axs1[2].autoscale()\n",
"fig1.tight_layout()\n",
"plt.show()\n",
"plt.pause(1)\n",
"\n", "\n",
"# loop through time steps of the pipeline\n", "# loop through time steps of the pipeline\n",
"for it_pipe in range(1,pipe.nt):\n", "for it_pipe in range(1,pipe.nt):\n",
"\n", "\n",
"# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n", "# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n",
" # set initial conditions for the reservoir time evolution calculted with e-RK4\n",
" V.pressure = p_old[0]\n", " V.pressure = p_old[0]\n",
" V.outflux = v_old[0]\n", " V.outflux = v_old[0]\n",
" # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n",
" for it_res in range(nt_eRK4):\n", " for it_res in range(nt_eRK4):\n",
" V.e_RK_4()\n", " V.e_RK_4() # call e-RK4 to update outflux\n",
" V.level = V.update_level(V.timestep)\n", " V.level = V.update_level(V.timestep) # \n",
" V.set_volume()\n", " V.set_volume() # update volume in reservoir\n",
" level_vec_2[it_res] = V.level\n", " level_vec_2[it_res] = V.level # save for plotting\n",
" if (V.level < critical_level_low) or (V.level > critical_level_high):\n", " if (V.level < critical_level_low) or (V.level > critical_level_high): # make sure to never exceed critical levels\n",
" i_max = it_pipe\n", " i_max = it_pipe # for plotting only calculated values\n",
" print('broke')\n", " break \n",
" break\n", " level_vec[it_pipe] = V.level \n",
" level_vec[it_pipe] = V.level\n",
"\n", "\n",
" # set boundary conditions for the next timestep of the characteristic method\n",
" p_boundary_res[it_pipe] = rho*g*V.level-v_old[1]**2*rho/2\n", " p_boundary_res[it_pipe] = rho*g*V.level-v_old[1]**2*rho/2\n",
" 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", " 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", " +dt*g*np.sin(alpha)\n",
"\n", "\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", " 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", " p_boundary_tur[it_pipe] = pipe.p_boundary_tur\n",
"\n", "\n",
" # perform the next timestep via the characteristic method\n",
" pipe.timestep_characteristic_method()\n", " pipe.timestep_characteristic_method()\n",
"\n", "\n",
"\n", " # plot some stuff\n",
" # lo_00.remove()\n", " # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n",
" # lo_01.remove()\n", " lo_00.remove()\n",
" lo_01.remove()\n",
" # lo_02.remove()\n", " # lo_02.remove()\n",
" # lo_00, = axs2[0].plot(pl_vec,pressure_conversion(pipe.p_old,'Pa','mWS')[0],marker='.',c='blue')\n", " # plot new pressure and velocity distribution in the pipeline\n",
" # lo_01, = axs2[1].plot(pl_vec,pipe.v_old,marker='.',c='blue')\n", " lo_00, = axs1[0].plot(pl_vec,pressure_conversion(pipe.p_old,'Pa','mWS')[0],marker='.',c='blue')\n",
" # lo_02, = axs2[2].plot(level_vec_2,c='blue')\n", " lo_01, = axs1[1].plot(pl_vec,pipe.v_old,marker='.',c='blue')\n",
" # fig2.suptitle(str(it_pipe))\n", " # lo_02, = axs1[2].plot(level_vec_2,c='blue')\n",
" # fig2.canvas.draw()\n", " fig1.suptitle(str(it_pipe))\n",
" # fig2.canvas.flush_events()\n", " fig1.canvas.draw()\n",
" # fig2.tight_layout()\n", " fig1.tight_layout()\n",
" # plt.pause(0.1) \n", " plt.pause(0.00001) \n",
"\n", "\n",
" # prepare for next loop\n",
" p_old = pipe.p_old\n", " p_old = pipe.p_old\n",
" v_old = pipe.v_old \n", " v_old = pipe.v_old \n",
"\n", "\n",
@@ -218,26 +279,31 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"%matplotlib qt5\n", "# plot time evolution of boundary pressure and velocity as well as the reservoir level\n",
"fig1,axs1 = plt.subplots(3,2)\n", "\n",
"axs1[0,0].plot(t_vec,pressure_conversion(p_boundary_res,'Pa','mWS')[0])\n", "fig2,axs2 = plt.subplots(3,2)\n",
"axs1[0,1].plot(t_vec,v_boundary_res)\n", "axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,'Pa','mWS')[0])\n",
"axs1[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,'Pa','mWS')[0])\n", "axs2[0,1].plot(t_vec,v_boundary_res)\n",
"axs1[1,1].plot(t_vec,v_boundary_tur)\n", "axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,'Pa','mWS')[0])\n",
"axs1[2,0].plot(t_vec,level_vec)\n", "axs2[1,1].plot(t_vec,v_boundary_tur)\n",
"axs1[0,0].set_title('Pressure Reservoir')\n", "axs2[2,0].plot(t_vec,level_vec)\n",
"axs1[0,1].set_title('Velocity Reservoir')\n", "axs2[0,0].set_title('Pressure reservoir')\n",
"axs1[1,0].set_title('Pressure Turbine')\n", "axs2[0,1].set_title('Velocity reservoir')\n",
"axs1[1,1].set_title('Velocity Turbine')\n", "axs2[1,0].set_title('Pressure turbine')\n",
"axs1[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2[1,1].set_title('Velocity turbine')\n",
"axs1[0,0].set_ylabel(r'$p$ [mWS]')\n", "axs2[2,0].set_title('Level reservoir')\n",
"axs1[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs1[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", "axs2[0,0].set_ylabel(r'$p$ [mWS]')\n",
"axs1[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs1[1,0].set_ylabel(r'$p$ [mWS]')\n", "axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
"axs1[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", "axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs1[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", "axs2[1,0].set_ylabel(r'$p$ [mWS]')\n",
"fig1.tight_layout()\n", "axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
"axs2[2,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
"axs2[2,0].set_ylabel(r'$h$ [m]')\n",
"axs2[2,1].axis('off')\n",
"fig2.tight_layout()\n",
"plt.show()" "plt.show()"
] ]
} }