end of day commit 27.06.2022

This commit is contained in:
Brantegger Georg
2022-06-27 16:16:24 +02:00
parent b6c01bfce2
commit 16bc55d47a
4 changed files with 169 additions and 59 deletions

View File

@@ -69,7 +69,6 @@ class Ausgleichsbecken_class:
def e_RK_4(self): def e_RK_4(self):
# Update to deal with non constant pipeline pressure!
yn = self.outflux/self.area_outflux yn = self.outflux/self.area_outflux
h = self.level h = self.level
dt = self.timestep dt = self.timestep

View File

@@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 34, "execution_count": 1,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -14,7 +14,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 35, "execution_count": 2,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -39,13 +39,13 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 36, "execution_count": 3,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
"data": { "data": {
"application/vnd.jupyter.widget-view+json": { "application/vnd.jupyter.widget-view+json": {
"model_id": "ece5839afa864ae2b836269b18279459", "model_id": "ad0ac3d52f5b4c5ba1d465ae850c9e69",
"version_major": 2, "version_major": 2,
"version_minor": 0 "version_minor": 0
}, },

View File

@@ -2,18 +2,19 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 37, "execution_count": 1,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"#imports\n", "#imports\n",
"import numpy as np\n", "import numpy as np\n",
"import matplotlib.pyplot as plt" "import matplotlib.pyplot as plt\n",
"from pressure_conversion import pressure_conversion"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 38, "execution_count": 2,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -21,21 +22,21 @@
"\n", "\n",
"g = 9.81 # gravitational acceleration [m/s²]\n", "g = 9.81 # gravitational acceleration [m/s²]\n",
"\n", "\n",
"L = 100 # length of pipeline [m]\n", "L = 1000 # length of pipeline [m]\n",
"rho = 1000 # density of water [kg/m³]\n", "rho = 1000 # density of water [kg/m³]\n",
"D = 1 # pipe diameter \n", "D = 1 # pipe diameter [m]\n",
"Q0 = 2 # initial flow in whole pipe [m³/s]\n", "Q0 = 2 # initial flow in whole pipe [m³/s]\n",
"h = 20 # water level in upstream reservoir [m]\n", "h = 20 # water level in upstream reservoir [m]\n",
"n = 10 # number of pipe segments in discretization\n", "n = 10 # number of pipe segments in discretization\n",
"nt = 150 # number of time steps\n", "nt = 11 # number of time steps\n",
"f_coeff = 0.1 # lambda = 0.01 Friction loss coefficient [m]\n", "f_D = 0.01 # Darcy friction factor\n",
"c = 400 # propagation velocity of the pressure wave\n", "c = 400 # propagation velocity of the pressure wave [m/s]\n",
"\n" "\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 39, "execution_count": 3,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -52,88 +53,121 @@
"\n", "\n",
"# storage vectors for old parameters\n", "# storage vectors for old parameters\n",
"v_old = np.full(nn,v0)\n", "v_old = np.full(nn,v0)\n",
"p_old = p0-(f_coeff*pl_vec/D*rho/2*v0**2) # ref Wikipedia: Rohrreibungszahls\n", "p_old = p0-(f_D*pl_vec/D*rho/2*v0**2) # ref Wikipedia: Rohrreibungszahls\n",
"\n", "\n",
"# storage vectors for new parameters\n", "# storage vectors for new parameters\n",
"v_new = np.zeros_like(v_old)\n", "v_new = np.zeros_like(v_old)\n",
"p_new = np.zeros_like(p_old)\n", "p_new = np.zeros_like(p_old)\n",
"\n", "\n",
"# storage vector for time evolution of parameters at node nn (at reservoir)\n", "# storage vector for time evolution of parameters at node 1 (at reservoir)\n",
"p_nn = np.zeros_like(t_vec)\n", "p_1 = np.zeros_like(t_vec)\n",
"v_nn = np.zeros_like(t_vec)\n", "v_1 = np.zeros_like(t_vec)\n",
"\n",
"# storage vector for time evolution of parameters at node N+1 (at valve)\n",
"p_np1 = np.zeros_like(t_vec)\n",
"v_np1 = np.zeros_like(t_vec)\n",
"\n" "\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 40, "execution_count": 4,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [
{
"data": {
"text/plain": [
"(-5.092958178940651, 5.092958178940651)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [ "source": [
"%matplotlib qt\n", "%matplotlib qt\n",
"# time loop\n", "# plotting preparation\n",
"\n", "\n",
"fig = plt.figure()\n", "fig1,axs1 = plt.subplots(2,1)\n",
"ax1 = fig.add_subplot(111)\n", "axs1[0].set_title('Pressure distribution in pipeline')\n",
"lo1, = ax1.plot(pl_vec,np.full_like(pl_vec,p0),marker='.')\n", "axs1[1].set_title('Velocity distribution in pipeline')\n",
"ax1.set_ylim([-20*p0,20*p0])\n",
"\n", "\n",
"lo_00, = axs1[0].plot(pl_vec,p_old,marker='.')\n",
"lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n",
"\n",
"axs1[0].set_ylim([-20*p0,20*p0])\n",
"axs1[1].set_ylim([-2*v0,2*v0])"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"8\n",
"[2.54647909 2.54647909 0.03242134 0.02836835 0.02431541 0.02026254\n",
" 0.01620977 0.01215711 0.00810458 0.0040522 0. ]\n",
"9\n",
"[2.54647909 0.03647353 0.03242052 0.02836756 0.02431467 0.02026188\n",
" 0.01620919 0.01215664 0.00810425 0.00405203 0. ]\n",
"10\n",
"[-2.46542799 -2.46568104 -2.46593345 -2.46618518 -2.4664362 -2.46668647\n",
" -2.46693595 -2.46718459 -2.46743236 -2.46767923 0. ]\n"
]
}
],
"source": [
"for it in range(nt):\n", "for it in range(nt):\n",
" # set boundary conditions\n", " # set boundary conditions\n",
" v_new[-1] = 0 # in front of the instantaneously closing valve, the velocity is 0\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", " p_new[0] = p0 # hydrostatic pressure from the reservoir\n",
"\n", "\n",
" # calculate the new parameters at first and last node\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", " v_new[0] = v_old[1]+1/(rho*c)*(p0-p_old[1])-f_D*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", " p_new[-1] = p_old[-2]+rho*c*v_old[-2]-rho*c*f_D*dt/(2*D) *abs(v_old[-2])*v_old[-2]\n",
"\n", "\n",
" # calculate parameters at second to second-to-last nodes \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", " #equation 2-30 plus 2-31 (and refactor for v_i^j+1) in block 2\n",
"\n", "\n",
" for i in range(1,nn-1):\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", " 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", " -f_D*dt/(4*D)*(abs(v_old[i-1])*v_old[i-1]+abs(v_old[i+1])*v_old[i+1])\n",
"\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", " 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", " -rho*c*f_D*dt/(4*D)*(abs(v_old[i-1])*v_old[i-1]-abs(v_old[i+1])*v_old[i+1])\n",
" \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", "\n",
" # store parameters of node 0 (at reservoir)\n", " lo_00.set_ydata(p_new)\n",
" p_nn[it] = p_old[-1]\n", " lo_01.set_ydata(v_new)\n",
" v_nn[it] = v_old[-1]\n", " \n",
" fig1.suptitle(str(it))\n",
" fig1.canvas.draw()\n",
" fig1.tight_layout()\n",
" plt.pause(0.2)\n",
"\n",
" # store parameters of node 1 (at reservoir)\n",
" p_1[it] = p_old[0]\n",
" v_1[it] = v_old[0]\n",
" # store parameters of node N+1 (at reservoir)\n",
" p_np1[it] = p_old[-1]\n",
" v_np1[it] = v_old[-1]\n",
"\n", "\n",
" # prepare for next loop\n", " # prepare for next loop\n",
" p_old = p_new\n", " p_old = p_new\n",
" v_old = v_new\n", " v_old = v_new\n",
" if it > 7:\n",
" print(it)\n",
" #print(pressure_conversion(p_new, input_unit= 'Pa', target_unit='Bar'))\n",
" print(v_new)\n",
"\n", "\n",
"\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": { "metadata": {

View File

@@ -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)