diff --git a/Druckrohrleitung/Druckrohrleitung_ETH_class_file.py b/Druckrohrleitung/Druckrohrleitung_ETH_class_file.py new file mode 100644 index 0000000..357186f --- /dev/null +++ b/Druckrohrleitung/Druckrohrleitung_ETH_class_file.py @@ -0,0 +1,171 @@ +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 + acceleration_unit = r'$\mathrm{m}/\mathrm{s}^2$' + angle_unit = '°' + area_unit = r'$\mathrm{m}^2$' + density_unit = r'$\mathrm{kg}/\mathrm{m}^3$' + flux_unit = r'$\mathrm{m}^3/\mathrm{s}$' + length_unit = 'm' + pressure_unit = 'Pa' + time_unit = 's' + velocity_unit = r'$\mathrm{m}/\mathrm{s}$' # for flux and pressure propagation + volume_unit = r'$\mathrm{m}^3$' + + acceleration_unit_print = 'm/s²' + angle_unit_print = '°' + area_unit_print = 'm²' + density_unit_print = 'kg/m³' + flux_unit_print = 'm³/s' + length_unit_print = 'm' + pressure_unit_print = 'Pa' + time_unit_print = 's' + velocity_unit_print = 'm/s' # for flux and pressure propagation + volume_unit_print = 'm³' + + +# init + def __init__(self,total_length,diameter,number_segments,pipeline_angle,Darcy_friction_factor,rho=1000,g=9.81): + self.length = total_length + self.dia = diameter + self.n_seg = number_segments + self.angle = pipeline_angle + self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient + self.density = 1000 + self.g = g + + self.dx = total_length/number_segments + self.l_vec = np.arange(0,(number_segments+1)*self.dx,self.dx) + + # initialize for get_info method + self.c = '--' + self.dt = '--' + +# setter + def set_pressure_propagation_velocity(self,c): + self.c = c + self.dt = self.dx/c + + def set_number_of_timesteps(self,number_timesteps): + self.nt = number_timesteps + if self.c == '--': + raise Exception('Please set the pressure propagation velocity before setting the number of timesteps.') + else: + self.t_vec = np.arange(0,self.nt*self.dt,self.dt) + + def set_initial_pressure(self,pressure,input_unit = 'Pa'): + p,_ = pressure_conversion(pressure,input_unit,target_unit=self.pressure_unit) + if np.size(p) == 1: + self.p0 = np.full_like(self.l_vec,p) + elif np.size(p) == np.size(self.l_vec): + self.p0 = p + else: + raise Exception('Unable to assign initial pressure. Input has to be of size 1 or' + np.size(self.l_vec)) + + #initialize the vectors in which the old and new pressures are stored for the method of characteristics + self.p_old = self.p0.copy() + self.p = np.empty_like(self.p_old) + + def set_initial_flow_velocity(self,velocity): + if np.size(velocity) == 1: + self.v0 = np.full_like(self.l_vec,velocity) + elif np.size(velocity) == np.size(self.l_vec): + self.v0 = velocity + else: + raise Exception('Unable to assign initial velocity. Input has to be of size 1 or' + np.size(self.l_vec)) + + #initialize the vectors in which the old and new velocities are stored for the method of characteristics + self.v_old = self.v0.copy() + self.v = np.empty_like(self.v_old) + + def set_boundary_conditions_next_timestep(self,v_reservoir,p_reservoir,v_turbine,input_unit_pressure = 'Pa'): + rho = self.density + c = self.c + f_D = self.f_D + dt = self.dt + D = self.dia + p_old = self.p_old[-2] # @ second to last node (the one before the turbine) + v_old = self.v_old[-2] # @ second to last node (the one before the turbine) + self.v_boundary_res = v_reservoir + self.v_boundary_tur = v_turbine + self.p_boundary_res,_ = pressure_conversion(p_reservoir,input_unit_pressure,target_unit=self.pressure_unit) + self.p_boundary_tur = p_old+rho*c*v_old-rho*c*f_D*dt/(2*D)*abs(v_old)*v_old + self.v[0] = self.v_boundary_res.copy() + self.v[-1] = self.v_boundary_tur.copy() + self.p[0] = self.p_boundary_res.copy() + self.p[-1] = self.p_boundary_tur.copy() + +# getter + def get_info(self): + new_line = '\n' + + # :<10 pads the self.value to be 10 characters wide + print_str = (f"The pipeline has the following attributes: {new_line}" + f"----------------------------- {new_line}" + f"Length = {self.length:<10} {self.length_unit_print} {new_line}" + f"Diameter = {self.dia:<10} {self.length_unit_print} {new_line}" + f"Number of segemnts = {self.n_seg:<10} {new_line}" + f"Number of nodes = {self.n_seg+1:<10} {new_line}" + f"Length per segment = {self.dx:<10} {self.length_unit_print} {new_line}" + f"Pipeline angle = {self.angle:<10} {self.angle_unit_print} {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"Pressure wave vel. = {self.c:<10} {self.velocity_unit_print} {new_line}" + f"Simulation timesteps = {self.dt:<10} {self.time_unit_print } {new_line}" + f"Number of timesteps = {self.nt:<10} {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") + + print(print_str) + + + def get_boundary_conditions_next_timestep(self,target_unit_pressure ='bar'): + print('The pressure at the reservoir for the next timestep is', '\n', \ + pressure_conversion(self.p_boundary_res,self.pressure_unit_print,target_unit_pressure), '\n', \ + 'The velocity at the reservoir for the next timestep is', '\n', \ + self.v_boundary_res, self.velocity_unit, '\n', \ + 'The pressure at the turbine for the next timestep is', '\n', \ + pressure_conversion(self.p_boundary_tur,self.pressure_unit_print,target_unit_pressure), '\n', \ + 'The velocity at the turbine for the next timestep is', '\n', \ + self.v_boundary_tur, self.velocity_unit) + + + def timestep_characteristic_method(self): + #number of nodes + nn = self.n_seg+1 + rho = self.density + c = self.c + f_D = self.f_D + dt = self.dt + D = self.dia + + 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]) \ + -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]) + + self.p[i] = 0.5*rho*c*(self.v_old[i-1]-self.v_old[i+1])+0.5*(self.p_old[i-1]+self.p_old[i+1]) \ + -rho*c*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]) + + self.p_old = self.p.copy() + self.v_old = self.v.copy() + + + + + + + + + + + diff --git a/Druckrohrleitung/Druckrohrleitung_class_file.py b/Druckrohrleitung/Druckrohrleitung_class_file.py index 357186f..815c49c 100644 --- a/Druckrohrleitung/Druckrohrleitung_class_file.py +++ b/Druckrohrleitung/Druckrohrleitung_class_file.py @@ -12,7 +12,7 @@ from functions.pressure_conversion import pressure_conversion class Druckrohrleitung_class: # units acceleration_unit = r'$\mathrm{m}/\mathrm{s}^2$' - angle_unit = '°' + angle_unit = 'rad' area_unit = r'$\mathrm{m}^2$' density_unit = r'$\mathrm{kg}/\mathrm{m}^3$' flux_unit = r'$\mathrm{m}^3/\mathrm{s}$' @@ -23,7 +23,7 @@ class Druckrohrleitung_class: volume_unit = r'$\mathrm{m}^3$' acceleration_unit_print = 'm/s²' - angle_unit_print = '°' + angle_unit_print = 'rad' area_unit_print = 'm²' density_unit_print = 'kg/m³' flux_unit_print = 'm³/s' @@ -89,21 +89,23 @@ class Druckrohrleitung_class: self.v = np.empty_like(self.v_old) def set_boundary_conditions_next_timestep(self,v_reservoir,p_reservoir,v_turbine,input_unit_pressure = 'Pa'): - rho = self.density - c = self.c - f_D = self.f_D - dt = self.dt - D = self.dia - p_old = self.p_old[-2] # @ second to last node (the one before the turbine) - v_old = self.v_old[-2] # @ second to last node (the one before the turbine) - self.v_boundary_res = v_reservoir - self.v_boundary_tur = v_turbine - self.p_boundary_res,_ = pressure_conversion(p_reservoir,input_unit_pressure,target_unit=self.pressure_unit) - self.p_boundary_tur = p_old+rho*c*v_old-rho*c*f_D*dt/(2*D)*abs(v_old)*v_old - self.v[0] = self.v_boundary_res.copy() - self.v[-1] = self.v_boundary_tur.copy() - self.p[0] = self.p_boundary_res.copy() - self.p[-1] = self.p_boundary_tur.copy() + rho = self.density + c = self.c + f_D = self.f_D + dt = self.dt + D = self.dia + g = self.g + alpha = self.angle + p_old = self.p_old[-2] # @ second to last node (the one before the turbine) + v_old = self.v_old[-2] # @ second to last node (the one before the turbine) + self.v_boundary_res = v_reservoir # at new timestep + self.v_boundary_tur = v_turbine # at new timestep + self.p_boundary_res,_ = pressure_conversion(p_reservoir,input_unit_pressure,target_unit=self.pressure_unit) + self.p_boundary_tur = p_old-rho*c*(v_turbine-v_old)+rho*c*dt*g*np.sin(alpha)-f_D*rho*c*dt/(2*D)*abs(v_old)*v_old + self.v[0] = self.v_boundary_res.copy() + self.v[-1] = self.v_boundary_tur.copy() + self.p[0] = self.p_boundary_res.copy() + self.p[-1] = self.p_boundary_tur.copy() # getter def get_info(self): @@ -142,19 +144,21 @@ class Druckrohrleitung_class: def timestep_characteristic_method(self): #number of nodes - nn = self.n_seg+1 - rho = self.density - c = self.c - f_D = self.f_D - dt = self.dt - D = self.dia + nn = self.n_seg+1 + rho = self.density + c = self.c + f_D = self.f_D + dt = self.dt + D = self.dia + g = self.g + alpha = self.angle 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]) \ - -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]) + 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]) - self.p[i] = 0.5*rho*c*(self.v_old[i-1]-self.v_old[i+1])+0.5*(self.p_old[i-1]+self.p_old[i+1]) \ - -rho*c*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]) + self.p[i] = 0.5*(self.p_old[i+1]+self.p_old[i-1]) - 0.5*rho*c*(self.v_old[i+1]-self.v_old[i-1]) \ + +f_D*rho*c*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]) self.p_old = self.p.copy() self.v_old = self.v.copy() diff --git a/Druckrohrleitung/Main_Programm.ipynb b/Druckrohrleitung/Main_Programm.ipynb index 83bafe6..fbcc88c 100644 --- a/Druckrohrleitung/Main_Programm.ipynb +++ b/Druckrohrleitung/Main_Programm.ipynb @@ -2,11 +2,12 @@ "cells": [ { "cell_type": "code", - "execution_count": 5, + "execution_count": 52, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", + "from numpy import sin, arcsin\n", "from Druckrohrleitung_class_file import Druckrohrleitung_class\n", "import matplotlib.pyplot as plt\n", "\n", @@ -21,7 +22,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 53, "metadata": {}, "outputs": [], "source": [ @@ -29,53 +30,56 @@ "#define constants\n", "\n", "g = 9.81 # gravitational acceleration [m/s²]\n", + "rho = 1000 # density of water [kg/m³]\n", "\n", "L = 1000 # length of pipeline [m]\n", - "rho = 1000 # density of water [kg/m³]\n", "D = 1 # pipe diameter [m]\n", "Q0 = 2 # initial flow in whole pipe [m³/s]\n", - "h = 20 # water level in upstream reservoir [m]\n", + "h_res = 20 # water level in upstream reservoir [m]\n", "n = 10 # number of pipe segments in discretization\n", - "nt = 500 # number of time steps after initial conditions\n", + "nt = 100 # 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 = 1e-5 # hydraulic head without reservoir [m] \n", + "alpha = arcsin(h_pipe/L) # Höhenwinkel der Druckrohrleitung \n", "\n", "\n", "# preparing the discretization and initial conditions\n", "\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", - "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*dt,dt) # time vector\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", + "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*dt,dt) # time vector\n", + "h_vec = np.arange(0,h_pipe+h_pipe/n,h_pipe/n) # hydraulic head of pipeline at each node\n", "\n", - "v0 = Q0/(D**2/4*np.pi)\n", - "p0 = (rho*g*h-v0**2*rho/2)\n", + "v_init = np.full(nn,Q0/(D**2/4*np.pi))\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", "# storage vectors for old parameters\n", - "v_old = np.full(nn,v0)\n", - "p_old = p0-(f_D*pl_vec/D*rho/2*v0**2) # ref Wikipedia: Darcy Weisbach\n", + "v_old = v_init.copy()\n", + "p_old = p_init.copy() \n", "\n", "# storage vectors for new parameters\n", - "v_new = np.zeros_like(v_old)\n", - "p_new = np.zeros_like(p_old)\n", + "v_new = np.empty_like(v_old)\n", + "p_new = np.empty_like(p_old)\n", "\n", - "# storage vector for time evolution of parameters at node 1 (at reservoir)\n", - "p_1 = np.full_like(t_vec,p0)\n", - "v_1 = np.full_like(t_vec,v0)\n", + "# storage vector for time evolution of parameters at node 0 (at reservoir)\n", + "p_0 = np.full_like(t_vec,p_init[0])\n", + "v_0 = np.full_like(t_vec,v_init[0])\n", "\n", "# storage vector for time evolution of parameters at node N+1 (at valve)\n", - "p_np1 = np.full_like(t_vec,p0)\n", - "v_np1 = np.full_like(t_vec,v0)\n", + "p_np1 = np.full_like(t_vec,p_init[-1])\n", + "v_np1 = np.full_like(t_vec,v_init[-1])\n", "\n", "for it in range(1,nt):\n", "\n", " # set boundary conditions\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", + " v_new[-1] = 0 # in front of the instantaneously closing valve, the velocity is 0\n", + " p_new[0] = p_init[0] # hydrostatic pressure from the reservoir\n", "\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_D*dt/(2*D)*abs(v_old[1])*v_old[1]\n", + " v_new[0] = v_old[1]+1/(rho*c)*(p_init[0]-p_old[1])+dt*g*sin(alpha)-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_D*dt/(2*D) *abs(v_old[-2])*v_old[-2]\n", "\n", " # calculate parameters at second to second-to-last nodes \n", @@ -83,7 +87,7 @@ "\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", - " -f_D*dt/(4*D)*(abs(v_old[i-1])*v_old[i-1]+abs(v_old[i+1])*v_old[i+1])\n", + " +dt*g*sin(alpha)-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", " 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_D*dt/(4*D)*(abs(v_old[i-1])*v_old[i-1]-abs(v_old[i+1])*v_old[i+1])\n", @@ -96,8 +100,8 @@ " v_old = v_new.copy()\n", "\n", " # store parameters of node 1 (at reservoir)\n", - " p_1[it] = p_new[0]\n", - " v_1[it] = v_new[0]\n", + " p_0[it] = p_new[0]\n", + " v_0[it] = v_new[0]\n", " # store parameters of node N+1 (at reservoir)\n", " p_np1[it] = p_new[-1]\n", " v_np1[it] = v_new[-1]" @@ -105,35 +109,7 @@ }, { "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "fig1,axs1 = plt.subplots(2,2)\n", - "axs1[0,0].plot(t_vec,p_1)\n", - "axs1[0,1].plot(t_vec,v_1)\n", - "axs1[1,0].plot(t_vec,p_np1)\n", - "axs1[1,1].plot(t_vec,v_np1)\n", - "axs1[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs1[0,0].set_ylabel(r'$p$ [Pa]')\n", - "axs1[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs1[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", - "axs1[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs1[1,0].set_ylabel(r'$p$ [Pa]')\n", - "axs1[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n", - "axs1[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n", - "\n", - "axs1[0,0].set_title('Pressure Reservoir')\n", - "axs1[0,1].set_title('Velocity Reservoir')\n", - "axs1[1,0].set_title('Pressure Turbine')\n", - "axs1[1,1].set_title('Velocity Turbine')\n", - "fig1.tight_layout()\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": 8, + "execution_count": 54, "metadata": {}, "outputs": [], "source": [ @@ -142,17 +118,17 @@ "pipe.set_pressure_propagation_velocity(c)\n", "pipe.set_number_of_timesteps(nt)\n", "\n", - "pipe.set_initial_pressure(p0)\n", - "pipe.set_initial_flow_velocity(v0)\n", - "pipe.set_boundary_conditions_next_timestep(v_1[0],p_1[0],v_np1[0])\n", + "pipe.set_initial_pressure(p_init)\n", + "pipe.set_initial_flow_velocity(v_init)\n", + "pipe.set_boundary_conditions_next_timestep(v_0[0],p_0[0],v_np1[0])\n", "\n", - "# storage vector for time evolution of parameters at node 1 (at reservoir)\n", - "pipe.p_1 = np.full_like(t_vec,p0)\n", - "pipe.v_1 = np.full_like(t_vec,v0)\n", + "# storage vector for time evolution of parameters at node 0 (at reservoir)\n", + "pipe.p_0 = np.full_like(t_vec,p_init[0])\n", + "pipe.v_0 = np.full_like(t_vec,v_init[0])\n", "\n", "# storage vector for time evolution of parameters at node N+1 (at valve)\n", - "pipe.p_np1 = np.full_like(t_vec,p0)\n", - "pipe.v_np1 = np.full_like(t_vec,v0)\n", + "pipe.p_np1 = np.full_like(t_vec,p_init[-1])\n", + "pipe.v_np1 = np.full_like(t_vec,v_init[-1])\n", "\n", "fig2,axs2 = plt.subplots(2,1)\n", "axs2[0].set_title('Pressure distribution in pipeline')\n", @@ -161,22 +137,22 @@ "axs2[0].set_ylabel(r'$p$ [Pa]')\n", "axs2[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs2[1].set_ylabel(r'$p$ [Pa]')\n", - "lo_00, = axs2[0].plot(pl_vec,pipe.p_old,marker='.')\n", + "lo_00, = axs2[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", - "axs2[0].set_ylim([-20*p0,20*p0])\n", - "axs2[1].set_ylim([-2*v0,2*v0])\n", + "axs2[0].set_ylim([-5*np.max(pressure_conversion(pipe.p_old,'Pa','mWs')[0]),5*np.max(pressure_conversion(pipe.p_old,'Pa','mWs')[0])])\n", + "axs2[1].set_ylim([-2*np.max(v_init),2*np.max(v_init)])\n", "fig2.tight_layout()\n", "\n", "\n", "for it in range(1,pipe.nt):\n", - " pipe.set_boundary_conditions_next_timestep(v_1[it],p_1[it],v_np1[it])\n", + " pipe.set_boundary_conditions_next_timestep(v_0[it],p_0[it],v_np1[it])\n", " pipe.timestep_characteristic_method()\n", " lo_00.set_ydata(pipe.p)\n", " lo_01.set_ydata(pipe.v)\n", "\n", - " # store parameters of node 1 (at reservoir)\n", - " pipe.p_1[it] = pipe.p[0]\n", - " pipe.v_1[it] = pipe.v[0]\n", + " # store parameters of node 0 (at reservoir)\n", + " pipe.p_0[it] = pipe.p[0]\n", + " pipe.v_0[it] = pipe.v[0]\n", " # store parameters of node N+1 (at reservoir)\n", " pipe.p_np1[it] = pipe.p[-1]\n", " pipe.v_np1[it] = pipe.v[-1]\n", @@ -184,18 +160,18 @@ " fig2.suptitle(str(it))\n", " fig2.canvas.draw()\n", " fig2.tight_layout()\n", - " plt.pause(0.001)\n" + " plt.pause(0.2)\n" ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 55, "metadata": {}, "outputs": [], "source": [ "fig3,axs3 = plt.subplots(2,2)\n", - "axs3[0,0].plot(t_vec,pipe.p_1)\n", - "axs3[0,1].plot(t_vec,pipe.v_1)\n", + "axs3[0,0].plot(t_vec,pipe.p_0)\n", + "axs3[0,1].plot(t_vec,pipe.v_0)\n", "axs3[1,0].plot(t_vec,pipe.p_np1)\n", "axs3[1,1].plot(t_vec,pipe.v_np1)\n", "axs3[0,0].set_title('Pressure Reservoir')\n", diff --git a/Messing Around/messy_nb.ipynb b/Messing Around/messy_nb.ipynb new file mode 100644 index 0000000..e69de29