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

fixed a coding mistake that lead to

Browse Source 
a missbehavior in the time evolution of the
reservoir
This commit is contained in:
Georg ´Brantegger
2022-07-28 16:26:04 +02:00
parent e90406d90b
commit dc5bcfe7f8
5 changed files with 135 additions and 2087 deletions
Download Patch File Download Diff File Expand all files Collapse all files

40
Ausgleichsbecken/Ausgleichsbecken_class_file.py
View File

@@ -73,10 +73,12 @@ class Ausgleichsbecken_class:
# sets the level in the reservoir and should only be called during initialization # sets the level in the reservoir and should only be called during initialization
if self.level == '--': if self.level == '--':
self.level = initial_level self.level = initial_level
self.volume = self.update_volume()
else: else:
raise Exception('Initial level was already set once. Use the .update_level(self,timestep) method to update level based on net flux.') raise Exception('Initial level was already set once. Use the .update_level(self,timestep) method to update level based on net flux.')
def set_level(self,level):
self.level = level
def set_influx(self,influx): def set_influx(self,influx):
# sets influx to the reservoir in m³/s # sets influx to the reservoir in m³/s
# positive influx means that liquid flows into the reservoir # positive influx means that liquid flows into the reservoir
@@ -95,15 +97,14 @@ class Ausgleichsbecken_class:
def set_pressure(self,pressure): def set_pressure(self,pressure):
# sets the static pressure present at the outlet of the reservoir # sets the static pressure present at the outlet of the reservoir
# units are used to convert and display the pressure
self.pressure = pressure self.pressure = pressure
def set_steady_state(self,ss_influx,ss_level,display_pressure_unit): def set_steady_state(self,ss_influx,ss_level,display_pressure_unit):
# set the steady state (ss) condition in which the net flux is zero # set the steady state (ss) condition in which the net flux is zero
# set pressure acting on the outflux area so that the level stays constant # set pressure acting on the outflux area so that the level stays constant
ss_outflux = ss_influx ss_outflux = ss_influx
ss_influx_vel = ss_influx/self.area ss_influx_vel = abs(ss_influx/self.area)
ss_outflux_vel = ss_outflux/self.area_outflux ss_outflux_vel = abs(ss_outflux/self.area_outflux)
ss_pressure = self.density*self.g*ss_level+self.density*ss_outflux_vel*(ss_influx_vel-ss_outflux_vel) ss_pressure = self.density*self.g*ss_level+self.density*ss_outflux_vel*(ss_influx_vel-ss_outflux_vel)
self.set_influx(ss_influx) self.set_influx(ss_influx)
@@ -115,6 +116,7 @@ class Ausgleichsbecken_class:
def get_info(self, full = False): def get_info(self, full = False):
new_line = '\n' new_line = '\n'
p = pressure_conversion(self.pressure,self.pressure_unit,self.pressure_unit_print) p = pressure_conversion(self.pressure,self.pressure_unit,self.pressure_unit_print)
outflux_vel = self.outflux/self.area_outflux
if full == True: if full == True:
@@ -129,7 +131,7 @@ class Ausgleichsbecken_class:
f"Volume in reservoir = {self.volume:<10} {self.volume_unit_print} {new_line}" f"Volume in reservoir = {self.volume:<10} {self.volume_unit_print} {new_line}"
f"Current influx = {self.influx:<10} {self.flux_unit_print} {new_line}" f"Current influx = {self.influx:<10} {self.flux_unit_print} {new_line}"
f"Current outflux = {self.outflux:<10} {self.flux_unit_print} {new_line}" f"Current outflux = {self.outflux:<10} {self.flux_unit_print} {new_line}"
f"Current outflux vel = {round(self.outflux_vel,3):<10} {self.velocity_unit_print} {new_line}" f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_print} {new_line}"
f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}" f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}"
f"Simulation timestep = {self.timestep:<10} {self.time_unit_print} {new_line}" f"Simulation timestep = {self.timestep:<10} {self.time_unit_print} {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}"
@@ -142,7 +144,7 @@ class Ausgleichsbecken_class:
f"Volume in reservoir = {self.volume:<10} {self.volume_unit_print} {new_line}" f"Volume in reservoir = {self.volume:<10} {self.volume_unit_print} {new_line}"
f"Current influx = {self.influx:<10} {self.flux_unit_print} {new_line}" f"Current influx = {self.influx:<10} {self.flux_unit_print} {new_line}"
f"Current outflux = {self.outflux:<10} {self.flux_unit_print} {new_line}" f"Current outflux = {self.outflux:<10} {self.flux_unit_print} {new_line}"
f"Current outflux vel = {round(self.outflux_vel,3):<10} {self.velocity_unit_print} {new_line}" f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_print} {new_line}"
f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}" f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_print} {new_line}"
f"----------------------------- {new_line}") f"----------------------------- {new_line}")
@@ -157,9 +159,6 @@ class Ausgleichsbecken_class:
def get_current_outflux(self): def get_current_outflux(self):
return self.outflux return self.outflux
def get_current_volume(self):
return self.volume
def get_current_pressure(self): def get_current_pressure(self):
return self.pressure return self.pressure
@@ -174,19 +173,14 @@ class Ausgleichsbecken_class:
# cannot set new level directly in this method, because it gets called to calcuate during the Runge Kutta # cannot set new level directly in this method, because it gets called to calcuate during the Runge Kutta
# to calculate a ficticious level at half the timestep # to calculate a ficticious level at half the timestep
net_flux = self.influx-self.outflux net_flux = self.influx-self.outflux
delta_V = net_flux*timestep delta_level = net_flux*timestep/self.area
new_level = (self.volume+delta_V)/self.area new_level = (self.level+delta_level)
return new_level return new_level
def update_volume(self):
# sets volume in reservoir based on self.level
return self.level*self.area
def update_pressure(self): def update_pressure(self):
influx_vel = self.influx/self.area influx_vel = abs(self.influx/self.area)
outflux_vel = self.outflux/self.area_outflux outflux_vel = abs(self.outflux/self.area_outflux)
p_new = self.density*self.g*self.level+self.density*outflux_vel*(influx_vel-outflux_vel) p_new = self.density*self.g*self.level+self.density*outflux_vel*(influx_vel-outflux_vel)
return p_new return p_new
def timestep_reservoir_evolution(self): def timestep_reservoir_evolution(self):
@@ -209,8 +203,10 @@ class Ausgleichsbecken_class:
ynp1 = yn + dt/6*(FODE_function(Y1,h,A,A_a,p,rho,g)+2*FODE_function(Y2,h_hs,A,A_a,p_hs,rho,g)+ \ ynp1 = yn + dt/6*(FODE_function(Y1,h,A,A_a,p,rho,g)+2*FODE_function(Y2,h_hs,A,A_a,p_hs,rho,g)+ \
2*FODE_function(Y3,h_hs,A,A_a,p_hs,rho,g)+ FODE_function(Y4,h,A,A_a,p,rho,g)) 2*FODE_function(Y3,h_hs,A,A_a,p_hs,rho,g)+ FODE_function(Y4,h,A,A_a,p,rho,g))
self.outflux = ynp1*A_a new_outflux = ynp1*A_a
self.level = self.update_level(dt) new_level = self.update_level(dt)
self.volume = self.update_volume() new_pressure = self.update_pressure()
self.pressure = self.update_pressure()
self.set_outflux(new_outflux)
self.set_level(new_level)
self.set_pressure(new_pressure)

16
Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 13, "execution_count": 1,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -21,7 +21,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 14, "execution_count": 2,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -51,7 +51,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 15, "execution_count": 3,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -71,14 +71,14 @@
"\n", "\n",
"# for while loop\n", "# for while loop\n",
"total_min_level = 0.01 # m\n", "total_min_level = 0.01 # m\n",
"total_max_time = 1 # s\n", "total_max_time = 100 # s\n",
"\n", "\n",
"nt = int(total_max_time//simulation_timestep)" "nt = int(total_max_time//simulation_timestep)"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 16, "execution_count": 4,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -119,7 +119,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 17, "execution_count": 5,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -149,7 +149,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 18, "execution_count": 6,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -158,7 +158,7 @@
"10.1" "10.1"
] ]
}, },
"execution_count": 18, "execution_count": 6,
"metadata": {}, "metadata": {},
"output_type": "execute_result" "output_type": "execute_result"
} }

4
Druckrohrleitung/Druckrohrleitung_class_file.py
View File

@@ -118,8 +118,8 @@ class Druckrohrleitung_class:
ss_v0 = np.full(self.n_seg+1,ss_flux/self.A) ss_v0 = np.full(self.n_seg+1,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 # the static pressure is given by static state pressure of the reservoir, corrected for the hydraulic head of the pipe and friction losses
ss_v_in_res = ss_flux/area_reservoir ss_v_in_res = abs(ss_flux/area_reservoir)
ss_v_out_res = ss_flux/self.A ss_v_out_res = abs(ss_flux/self.A)
ss_pressure_res = self.density*self.g*(ss_level_reservoir)+self.density*ss_v_out_res*(ss_v_in_res-ss_v_out_res) ss_pressure_res = self.density*self.g*(ss_level_reservoir)+self.density*ss_v_out_res*(ss_v_in_res-ss_v_out_res)
ss_pressure = ss_pressure_res+(self.density*self.g*h_vec)-(self.f_D*pl_vec/self.dia*self.density/2*ss_v0**2) ss_pressure = ss_pressure_res+(self.density*self.g*h_vec)-(self.f_D*pl_vec/self.dia*self.density/2*ss_v0**2)

2054
Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb
View File

File diff suppressed because it is too large Load Diff

100
Untertweng.ipynb
View File

@@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 1,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -17,7 +17,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 2,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -26,7 +26,7 @@
"#Turbine\n", "#Turbine\n",
"Q_nenn = 0.85 # m³/s\n", "Q_nenn = 0.85 # m³/s\n",
"p_nenn = pressure_conversion(10.6,'bar','Pa')\n", "p_nenn = pressure_conversion(10.6,'bar','Pa')\n",
"closing_time = 30 #s\n", "closing_time = 5 #s\n",
"\n", "\n",
"# physics\n", "# physics\n",
"g = 9.81 # gravitational acceleration [m/s²]\n", "g = 9.81 # gravitational acceleration [m/s²]\n",
@@ -38,14 +38,14 @@
"A_pipe = D**2/4*np.pi # pipeline area\n", "A_pipe = D**2/4*np.pi # pipeline area\n",
"h_pipe = 105 # 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", "alpha = np.arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n",
"n = 200 # number of pipe segments in discretization\n", "n = 50 # number of pipe segments in discretization\n",
"# consider replacing Q0 with a vector be be more flexible in initial conditions\n", "# consider replacing Q0 with a vector be be more flexible in initial conditions\n",
"# Q0 = Q_nenn # 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", "# v0 = Q0/A_pipe # initial flow velocity [m/s]\n",
"f_D = 0.014 # Darcy friction factor\n", "f_D = 0.014 # Darcy friction factor\n",
"c = 500. # 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", "# consider prescribing a total simulation time and deducting the number of timesteps from that\n",
"nt = 20000 # number of time steps after initial conditions\n", "nt = 2500 # number of time steps after initial conditions\n",
"\n", "\n",
"# derivatives of the pipeline constants\n", "# derivatives of the pipeline constants\n",
"dx = L/n # length of each pipe segment\n", "dx = L/n # length of each pipe segment\n",
@@ -69,7 +69,7 @@
"critical_level_high = np.inf # 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 = 100 # number of simulation steps of reservoir in between timesteps of pipeline \n",
"simulation_timestep = dt/nt_eRK4\n", "simulation_timestep = dt/nt_eRK4\n",
"\n", "\n",
"\n" "\n"
@@ -77,7 +77,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 3,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -90,7 +90,7 @@
"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_steady_state(initial_flux,initial_level,pl_vec,h_vec)\n", "pipe.set_steady_state(initial_flux,initial_level,area_base,pl_vec,h_vec)\n",
"\n", "\n",
"initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n", "initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n",
"\n", "\n",
@@ -105,7 +105,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 4,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -142,7 +142,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 5,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -170,21 +170,20 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 6,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# loop through time steps of the pipeline\n", "# loop through time steps of the pipeline\n",
"for it_pipe in range(1,pipe.nt+1):\n", "for it_pipe in range(1,pipe.nt+1):\n",
"\n", "\n",
" if t_vec[it_pipe]>20:\n", " if it_pipe == 250:\n",
" if V.get_current_influx() > 0:\n", " V.set_influx(0.)\n",
" V.set_influx(np.max([V.get_current_influx()-initial_flux*5*1e-3,0.]))\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", " # set initial conditions for the reservoir time evolution calculted with e-RK4\n",
" V.set_pressure = p_old[0]\n", " V.set_pressure(p_old[0])\n",
" V.set_outflux = v_old[0]*area_outflux\n", " V.set_outflux(v_old[0]*area_outflux)\n",
" # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\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.timestep_reservoir_evolution() \n", " V.timestep_reservoir_evolution() \n",
@@ -233,7 +232,24 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0\n"
]
}
],
"source": [
"print(V.get_current_influx())"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -274,11 +290,57 @@
"fig2.tight_layout()\n", "fig2.tight_layout()\n",
"plt.show()" "plt.show()"
] ]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The cuboid reservoir has the following attributes: \n",
"----------------------------- \n",
"Base area = 74.0 m² \n",
"Outflux area = 0.636 m² \n",
"Current level = 7.875725956447418 m\n",
"Critical level low = 0.0 m \n",
"Critical level high = inf m \n",
"Volume in reservoir = -- m³ \n",
"Current influx = 0.0 m³/s \n",
"Current outflux = -0.1415386124341686 m³/s \n",
"Current outflux vel = -0.222 m/s \n",
"Current pipe pressure = 0.772 bar \n",
"Simulation timestep = 0.0004052 s \n",
"Density of liquid = 1000 kg/m³ \n",
"----------------------------- \n",
"\n",
"9.22707730779877\n",
"10.57842865915012\n",
"11.92978001050147\n",
"13.281131361852822\n",
"14.632482713204173\n",
"15.983834064555523\n",
"17.335185415906874\n",
"18.686536767258225\n",
"20.037888118609576\n",
"21.389239469960927\n"
]
}
],
"source": [
"V.get_info(full=True)\n",
"V.set_outflux(-10.)\n",
"for i in range(10):\n",
" V.level = V.update_level(10.)\n",
" print(V.get_current_level())"
]
} }
], ],
"metadata": { "metadata": {
"kernelspec": { "kernelspec": {
"display_name": "Python 3.8.13 ('Georg_DT_Slot3')", "display_name": "Python 3.8.13 ('DT_Slot_3')",
"language": "python", "language": "python",
"name": "python3" "name": "python3"
}, },
@@ -297,7 +359,7 @@
"orig_nbformat": 4, "orig_nbformat": 4,
"vscode": { "vscode": {
"interpreter": { "interpreter": {
"hash": "84fb123bdc47ab647d3782661abcbe80fbb79236dd2f8adf4cef30e8755eb2cd" "hash": "4a28055eb8a3160fa4c7e4fca69770c4e0a1add985300856aa3fcf4ce32a2c48"
} }
} }
}, },