Preparation for non static pipeline pressure

static_pipeline_pressure/Ausgleichsbecken.py

@@ -0,0 +1,67 @@
+import numpy as np 
+
+def Volume_trend(influx, outflux, timestep=1, V_0=0):
+    '''
+    Returns the trend and the volume and the final volume,  defined
+    by influx and outflux patterns. The optional parameter timestep
+    defines the time increment over which the fluxes are changing.
+    '''
+    net_flux = influx-outflux
+    delta_V = net_flux*timestep
+    V_trend = V_0+np.cumsum(delta_V)
+    V_end = V_trend[-1]
+    return V_end,  V_trend
+
+def Height_trend(V_trend, area=1, h_crit_low=-np.inf, h_crit_high=np.inf):
+    '''
+    Returns the trend and the height and the final height,  defined
+    by influx and outflux patterns as well as the crosssection area.
+    The optional parameters h_crit_low/high indicate limits that the height
+    should never exceed. If this occures,  TRUE is returned in the corresponding
+    h_crit_flag.
+    '''
+    h_trend = V_trend/area
+    h_crit_flag_low = np.any(h_trend <= h_crit_low)
+    h_crit_flag_high = np.any(h_trend >= h_crit_high)
+    h_end = h_trend[-1]
+    return h_trend, h_end, h_crit_flag_low, h_crit_flag_high
+
+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, p)+2*FODE_function(Y2, h_hs, alpha, p_hs)+ \
+        2*FODE_function(Y3, h_hs, alpha, p_hs)+ FODE_function(Y4, h, alpha, p))
+
+
+
+
+## testing
+# if __name__ == "__main__":
+#     influx = np.full([1, 100],  6)
+#     outflux = np.full_like(influx,  4)
+#     V_end,  V_trend = Volume_trend(influx,  outflux, timestep=0.5, V_0 = 100)
+#     print(V_end)
+#     print(V_trend)

static_pipeline_pressure/Ausgleichsbecken_class_file.py

@@ -0,0 +1,87 @@
+from Ausgleichsbecken import FODE_function, get_h_halfstep, get_p_halfstep
+from functions.pressure_conversion import pressure_conversion
+class Ausgleichsbecken_class:
+# units
+    area_unit           = r'$\mathrm{m}^2$'
+    area_outflux_unit   = r'$\mathrm{m}^2$'
+    level_unit          = 'm'
+    volume_unit         = r'$\mathrm{m}^3$'
+    flux_unit           = r'$\mathrm{m}^3/\mathrm{s}$'
+    time_unit           = 's'
+    pressure_unit       = 'Pa'
+
+# init
+    def __init__(self,area,outflux_area,level_min,level_max,timestep = 1):
+        self.area           = area              # base area of the rectangular structure
+        self.area_outflux   = outflux_area      # area of the outlet towards the pipeline        
+        self.level_min      = level_min         # lowest  allowed water level    
+        self.level_max      = level_max         # highest allowed water level
+        self.timestep       = timestep          # timestep of the simulation    
+
+# setter
+    def set_volume(self):
+        self.volume = self.level*self.area 
+
+    def set_initial_level(self,initial_level):
+        self.level = initial_level
+        self.set_volume()
+
+    def set_influx(self,influx):
+        self.influx = influx
+
+    def set_outflux(self,outflux):
+        self.outflux = outflux
+
+# getter
+    def get_area(self):
+        print('The base area of the cuboid reservoir is', self.area, self.area_unit)
+
+    def get_outflux_area(self):
+        print('The outflux area from the cuboid reservoir to the pipeline is', \
+            self.area_outflux, self.area_outflux_unit)
+    
+    def get_level(self):
+        print('The current level in the reservoir is', self.level , self.level_unit)
+
+    def get_crit_levels(self):
+        print('The critical water levels in the reservoir are:  \n',\
+            '   Minimum:', self.level_min , self.level_unit , '\n',\
+            '   Maximum:', self.level_max , self.level_unit )
+
+    def get_volume(self):
+        print('The current water volume in the reservoir is', self.volume, self.volume_unit)
+
+    def get_timestep(self):
+        print('The timestep for the simulation is' , self.timestep, self.time_unit)
+
+    def get_influx(self):
+        print('The current influx is', self.influx, self.flux_unit)
+
+    def get_outflux(self):
+        print('The current outflux is', self.outflux, self.flux_unit)
+
+# methods
+    def update_level(self,timestep):
+        net_flux = self.influx-self.outflux
+        delta_V = net_flux*timestep
+        new_level = (self.volume+delta_V)/self.area
+        return new_level
+
+
+    def e_RK_4(self):
+        # Update to deal with non constant pipeline pressure!
+        yn = self.outflux/self.area_outflux
+        h = self.level
+        dt = self.timestep
+        p,_ = pressure_conversion(self.initial_pressure,self.pressure_unit,'Pa')
+        p_hs,_ = pressure_conversion(self.initial_pressure,self.pressure_unit,'Pa')
+        alpha = (self.area_outflux/self.area-1)
+        h_hs = self.update_level(dt/2)
+        Y1 = yn
+        Y2 = yn + dt/2*FODE_function(Y1, h, alpha, self.initial_pressure)
+        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, p)+2*FODE_function(Y2, h_hs, alpha, p_hs)+ \
+            2*FODE_function(Y3, h_hs, alpha, p_hs)+ FODE_function(Y4, h, alpha, p))
+
+        self.outflux = ynp1*self.area_outflux