From 50fb2d4a0429dc0851fac3c03232774aca458ddb Mon Sep 17 00:00:00 2001 From: Brantegger Georg Date: Mon, 20 Jun 2022 14:28:49 +0200 Subject: [PATCH] Preparation for non static pipeline pressure --- static_pipeline_pressure/Ausgleichsbecken.py | 67 ++++++++++++++ .../Ausgleichsbecken_class_file.py | 87 +++++++++++++++++++ .../e-RK4-Test.ipynb | 0 3 files changed, 154 insertions(+) create mode 100644 static_pipeline_pressure/Ausgleichsbecken.py create mode 100644 static_pipeline_pressure/Ausgleichsbecken_class_file.py rename e-RK4-Test.ipynb => static_pipeline_pressure/e-RK4-Test.ipynb (100%) diff --git a/static_pipeline_pressure/Ausgleichsbecken.py b/static_pipeline_pressure/Ausgleichsbecken.py new file mode 100644 index 0000000..e887617 --- /dev/null +++ b/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) diff --git a/static_pipeline_pressure/Ausgleichsbecken_class_file.py b/static_pipeline_pressure/Ausgleichsbecken_class_file.py new file mode 100644 index 0000000..a6e2500 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/e-RK4-Test.ipynb b/static_pipeline_pressure/e-RK4-Test.ipynb similarity index 100% rename from e-RK4-Test.ipynb rename to static_pipeline_pressure/e-RK4-Test.ipynb