import math import matplotlib as mpl import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation from collections import deque import datetime import time from pprint import pprint import csv ### Create Surface & Load Files ### import tkinter as tk import tkinter.filedialog as fd from time import sleep root = tk.Tk() root.withdraw() surface = [] taper_setup = {} taper_setup['tubing_head_pressure'] = 50 taper_setup['fluid_gradient'] = .45 taper_setup['stuffing_box_friction']= 100 taper_setup['c'] = [0, .08 , .08] taper_setup['rod_length_data'] = [0, 10095, 300 ] taper_setup['num_tapers'] = len(taper_setup['rod_length_data'])-1 taper_setup['rod_diameter_data'] = [0, .750 , 1.5 ] taper_setup['rod_youngs_modulus_data'] = [0, 30.5, 30.5 ] taper_setup['rod_weight_data'] = [0, 1.634 , 6] taper_setup['dt'] = 0.0500 taper_setup['tubing_id'] = 1.995 taper_setup['tubing_od'] = 2.375 unit_config = {} unit_config['anchor_depth'] = 8923.0 unit_config['well_name'] = "Python Demo Well" unit_config['pump_diameter'] = 1.25 unit_config['pump_area'] = (unit_config['pump_diameter']**2) * math.pi load_before = 0 load_after = 0 load_before_3 = 0 load_after_3 = 0 blank_array = [0]*100 top_pos_array = deque() top_load_array = deque() for i in range(0,10): top_pos_array.append(blank_array[:]) top_load_array.append(blank_array[:]) pump_pos_array = [] pump_load_array = [] pr_l = [0,0] pr_p = [0,0] pmp_p = [0,0] pmp_l = [0,0] count = [0]* (taper_setup['num_tapers']+1) position_load_counter = 0 surface_pos_array = [] surface_load_array = [] class Card: def __init__(self, id, s_p_array, s_l_array, d_p_array, d_l_array, t, u_c): self.card_id = id self.date_time = datetime.datetime.now() self.tapers = t self.unit_config = u_c self.surface_pos = s_p_array[:] self.surface_load = s_l_array[:] self.downhole_pos = d_p_array[:] self.downhole_load = d_l_array[:] self.downhole_fluid_load_adjustment = 223.907 print("downhole_fluid_load_adjustment not configured") self.card_type = "NotConfigured" print("card_type not configured") def set_strokes_per_minute(self, strokes_per_minute): self.strokes_per_minute = strokes_per_minute def __find_nearest_point_to_load(self, direction, card_half ,target_load, pos_array, load_array): left_ignore = min(pos_array) + (min(pos_array)+ max(pos_array))*(.75) right_ignore = min(pos_array) + (min(pos_array)+ max(pos_array))*(.25) if (direction=="up"): for i in range(1, len(pos_array)): if (card_half=="right" and pos_array[i] >= right_ignore): if ((load_array[i] <= target_load) & (load_array[i-1] > target_load)): return [pos_array[i], load_array[i]] elif( card_half =="left" and pos_array[i] <= left_ignore): if ((load_array[i] >= target_load) & (load_array[i-1] < target_load)): return [pos_array[i], load_array[i]] else: for j in range(len(pos_array)-1, 1, -1): if (card_half=="right" and pos_array[i] >= right_ignore): if ((load_array[i] >= target_load) & (load_array[i-1] < target_load)): return [pos_array[i], load_array[i]] elif( card_half =="left" and pos_array[i] <= left_ignore): if ((load_array[i] <= target_load) & (load_array[i-1] > target_load)): return [pos_array[i], load_array[i]] def __find_incremental_load(self, target_position, pos_array, load_array): for i in range(1, len(pos_array)): if ((pos_array[i] > target_position) and (pos_array[i-1] < target_position)): up_point_greater = [pos_array[i], load_array[i]] up_point_lesser = [pos_array[i-1], load_array[i-1]] if ((pos_array[i] < target_position) and (pos_array[i-1] > target_position)): down_point_greater = [pos_array[i], load_array[i]] down_point_lesser = [pos_array[i-1], load_array[i-1]] m_up = (up_point_greater[1] - up_point_lesser[1]) / (up_point_greater[0] - up_point_lesser[0]) b_up = up_point_greater[1] - (m_up * up_point_greater[0]) up_middle = (m_up*target_position) + b_up m_down = (down_point_greater[1] - down_point_lesser[1]) / (down_point_greater[0] - down_point_lesser[0]) b_down = down_point_greater[1] - (m_down * down_point_greater[0]) down_middle = (m_down*target_position) + b_down return up_middle - down_middle def find_limits(self): surface_max_position_position = max(self.surface_pos) surface_max_position_load = self.surface_load[self.surface_pos.index(surface_max_position_position)] self.surface_max_position = [surface_max_position_position, surface_max_position_load] surface_min_position_position = min(self.surface_pos) surface_min_position_load = self.surface_load[self.surface_pos.index(surface_min_position_position)] self.surface_min_position = [surface_min_position_position, surface_min_position_load] surface_max_load_load = max(self.surface_load) surface_max_load_position = self.surface_pos[self.surface_load.index(surface_max_load_load)] self.surface_max_load = [surface_max_load_position, surface_max_load_load] surface_min_load_load = min(self.surface_load) surface_min_load_position = self.surface_pos[self.surface_load.index(surface_min_load_load)] self.surface_min_load = [surface_min_load_position, surface_min_load_load] downhole_max_load_load = max(self.downhole_load) downhole_max_load_position = self.downhole_pos[self.downhole_load.index(downhole_max_load_load)] self.downhole_max_load = [downhole_max_load_position, downhole_max_load_load] downhole_min_load_load = min(self.downhole_load) downhole_min_load_position = self.downhole_pos[self.downhole_load.index(downhole_min_load_load)] self.downhole_min_load = [downhole_min_load_position, downhole_min_load_load] downhole_min_position_position = min(self.downhole_pos) downhole_min_position_load = self.downhole_load[self.downhole_pos.index(downhole_min_position_position)] self.downhole_min_position = [downhole_min_position_position, downhole_min_position_load] downhole_max_position_position = max(self.downhole_pos) downhole_max_position_load = self.downhole_load[self.downhole_pos.index(downhole_max_position_position)] self.downhole_max_position = [downhole_max_position_position, downhole_max_position_load] self.downhole_fluid_load = (self.downhole_max_load[1] - self.downhole_min_load[1]) - self.downhole_fluid_load_adjustment self.surface_stroke_length = self.surface_max_position[0] - self.surface_min_position[0] def print_limits(self): if not ('self.surface_stroke_length' in locals()): self.find_limits() print("Surface Max Position:", self.surface_max_position) print("Surface Min Position:", self.surface_min_position) print("Surface Max Load:", self.surface_max_load) print("Surface Min Load:", self.surface_min_load) print("Downhole Max Load:", self.downhole_max_load) print("Downhole Min Load:", self.downhole_min_load) print("Downhole Min Position:", self.downhole_min_position) print("Downhole Max Position:", self.downhole_max_position) def find_corners(self, pct_top, pct_bottom): if not ('self.surface_stroke_length' in locals()): self.find_limits() downhole_load_span = self.downhole_max_load[1] - self.downhole_min_load[1] half_load_load = self.downhole_min_load[1] + downhole_load_span / 2 half_load = self.__find_nearest_point_to_load("up", "left", half_load_load , self.downhole_pos, self.downhole_load) if(self.downhole_min_position[0] < half_load[0]): self.corner_bottom_left = self.downhole_min_position self.corner_top_left = self.__find_nearest_point_to_load("up", "left", self.downhole_max_load[1]-downhole_load_span*(pct_top/100), self.downhole_pos, self.downhole_load) else: self.corner_bottom_left = self.__find_nearest_point_to_load("up", "left", self.downhole_min_load[1]+downhole_load_span*(pct_bottom/100), self.downhole_pos, self.downhole_load) self.corner_top_left = self.downhole_min_position self.corner_top_right = self.downhole_max_position self.corner_fillage = self.__find_nearest_point_to_load("up", "right", self.downhole_min_load[1]+downhole_load_span*(pct_bottom/100), self.downhole_pos, self.downhole_load) self.corner_full_fillage_point = [self.corner_bottom_left[0]+ self.corner_top_right[0] - self.corner_top_left[0], self.corner_bottom_left[1]] def print_corners(self): if not ('self.corner_fillage' in locals()): self.find_corners(25, 25) print("Bottom Left:", self.corner_bottom_left) print("Bottom Left:", self.corner_top_left) print("Bottom Left:", self.corner_top_right) print("Bottom Left:", self.corner_fillage) print("Bottom Left:", self.corner_full_fillage_point) def calc_fillage(self, pct_top, pct_bottom): if not ('self.corner_fillage' in locals()): self.find_corners(pct_top, pct_bottom) self.downhole_gross_stroke = self.downhole_max_position[0] - self.downhole_min_position[0] self.tubing_movement = 12 * (tapers.rod_depth_total - self.unit_config['anchor_depth'])*self.downhole_fluid_load / (30500000* self.tapers.tubing_cross_sectional_area) self.downhole_adjusted_gross_stroke = self.downhole_gross_stroke - self.tubing_movement self.downhole_net_stroke = self.corner_fillage[0] - self.corner_bottom_left[0] self.fillage_estimated = ((self.corner_fillage[0] - self.corner_bottom_left[0])/(self.corner_full_fillage_point[0] - self.corner_bottom_left[0])) * 100 self.fillage_calculated = ((self.downhole_net_stroke + self.tubing_movement) / self.downhole_gross_stroke) * 100 def listSurface(self): for i in range(0,len(self.surface_pos)): print([self.surface_pos[i], self.surface_load[i]]) def listDownhole(self): for i in range(0,len(self.downhole_pos)): print([self.downhole_pos[i], self.downhole_load[i]]) def calc_hp(self, riemann_slices): self.polished_rod_horsepower = 0 self.pump_horsepower = 0 downSlices = [] if not ('self.surface_stroke_length' in locals()): self.find_limits() if not ('self.downhole_net_stroke' in locals()): self.calc_fillage(25,25) dx_surface = self.surface_stroke_length / (riemann_slices + 1) dx_downhole = self.downhole_net_stroke / (riemann_slices + 1) for i in range(1, riemann_slices): target_surface = dx_surface * i + self.surface_min_position[0] target_downhole = dx_downhole*i + self.downhole_min_position[0] slice_surface = self.__find_incremental_load(target_surface, self.surface_pos, self.surface_load) slice_downhole = self.__find_incremental_load(target_downhole, self.downhole_pos, self.downhole_load) self.polished_rod_horsepower += (dx_surface / 12) * slice_surface * self.strokes_per_minute / 33000 self.pump_horsepower += (dx_downhole / 12) * slice_downhole * self.strokes_per_minute / 33000 def calc_fluid_level(self): if not ('self.downhole_fluid_load' in locals()): self.find_limits() self.pump_intake_pressure = self.tapers.params['fluid_gradient'] * self.tapers.rod_depth_total - (self.downhole_fluid_load / self.unit_config['pump_area']) self.fluid_above_pump = self.pump_intake_pressure / self.tapers.params['fluid_gradient'] def write_csv(self): if not ('self.fluid_above_pump' in locals()): self.calc_fluid_level() if not ('self.polished_rod_horsepower' in locals()): self.calc_hp(100) if not ('self.fillage_calculated' in locals()): self.calc_fillage(25,25) csvData = {} csvData['card_id'] = self.card_id csvData['num_tapers'] = self.tapers.params['num_tapers'] csvData['num_points'] = len(self.surface_pos) csvData['card_type'] = self.card_type csvData['_year'] = self.date_time.year csvData['_month'] = self.date_time.month csvData['_day'] = self.date_time.day csvData['_hour'] = self.date_time.hour csvData['_minute'] = self.date_time.minute csvData['_second'] = self.date_time.second csvData['well_name'] = self.unit_config['well_name'] csvData['tubing_head_pressure'] = self.tapers.params['tubing_head_pressure'] csvData['fluid_gradient'] = self.tapers.params['fluid_gradient'] csvData['stuffing_box_friction'] = self.tapers.params['stuffing_box_friction'] csvData['dt'] = self.tapers.dt csvData['downhole_max_load'] = self.downhole_max_load[1] csvData['downhole_min_load'] = self.downhole_min_load[1] csvData['downhole_max_position'] = self.downhole_max_position[0] csvData['downhole_min_position'] = self.downhole_min_position[0] csvData['downhole_gross_stroke'] = self.downhole_gross_stroke csvData['downhole_adjusted_gross_stroke'] = self.downhole_adjusted_gross_stroke csvData['downhole_net_stroke'] = self.downhole_net_stroke csvData['downhole_fluid_load'] = self.downhole_fluid_load csvData['surface_max_load'] = self.surface_max_load[1] csvData['surface_min_load'] = self.surface_min_load[1] csvData['surface_max_position'] = self.surface_max_position[0] csvData['surface_min_position'] = self.surface_min_position[0] csvData['tubing_movement'] = self.tubing_movement csvData['surface_stroke_length'] = self.surface_stroke_length csvData['fillage_percent'] = self.fillage_calculated csvData['polished_rod_horsepower'] = self.polished_rod_horsepower csvData['pump_horsepower'] = self.pump_horsepower csvData['strokes_per_minute'] = self.strokes_per_minute csvData['fluid_above_pump'] = self.fluid_above_pump file_date = str(self.date_time.year) + str(self.date_time.month) + str(self.date_time.day) file_time = str(self.date_time.hour) + str(self.date_time.minute) + str(self.date_time.second) file_fillage = str(round(self.fillage_calculated,3)) filename = file_date + '_' + file_time + '_' + str(self.card_id) + '_' + self.card_type + '_' + str(self.fillage_calculated).replace(".", "-") + ".csv" with open(filename, 'wt', newline='') as card_file: writer = csv.writer(card_file) for data_point in csvData: writer.writerow([data_point, csvData[data_point]]) writer.writerow(["s_pos", "s_load"]) for i in range(0, len(self.surface_pos)): writer.writerow([self.surface_pos[i], self.surface_load[i]]) writer.writerow(["d_pos", "d_load"]) for j in range(0, len(self.downhole_pos)): writer.writerow([self.downhole_pos[j], self.downhole_load[j]]) def plot(self): fig = plt.figure(figsize=(12, 9)) fig.canvas.set_window_title('Well Load & Position Cards') ax1 = fig.add_subplot(2,1,1) ax2 = fig.add_subplot(2,1,2) ax1.set_title("Surface Card") ax2.set_title("Downhole Card") ax1.set_xlabel('Position') ax2.set_xlabel('Position') ax1.set_ylabel('Load') ax2.set_ylabel('Load') ax1.fill(self.surface_pos, self.surface_load, 'g') ax2.fill(self.downhole_pos, self.downhole_load, 'b') leg1 = ax1.legend() leg2 = ax2.legend() fig.subplots_adjust(hspace=.4) ax1.grid(True) ax2.grid(True) plt.ion() plt.show() class Taper: def __init__(self,params): self.params = params self.buoyant_force_total = 0 self.rod_weight_in_fluid_total = 0 self.rod_weight_in_air_total = 0 self.weight_data_total = 0 self.rod_depth_total = 0 self.annular_force_data_total = 0 self.a = [] self.area = [] self.pressure = [self.params['tubing_head_pressure']] self.buoyant_force = [] self.stretch = [] self.weight_data = [] self.annular_force_data = [] self.force = [] self.alpha = [] self.x_over_a = [] self.factor_array = [] self.lag_index_array = [] self.center_point = [] self.sum_center_point = [] self.length_required =[] self.dt = params['dt'] self.rod_depth = [] self.rod_weight_in_air = [] self.rod_weight_in_fluid = [] for x in range(0, self.params['num_tapers'] + 2): self.area.append(0) self.pressure.append(0) self.buoyant_force.append(0) self.rod_weight_in_air.append(0) self.rod_depth.append(0) self.rod_weight_in_fluid.append(0) self.weight_data.append(0) self.annular_force_data.append(0) self.a.append(0) self.stretch.append(0) self.force.append(0) self.alpha.append(0) self.x_over_a.append(0) self.factor_array.append(0) self.lag_index_array.append(0) self.center_point.append(0) self.sum_center_point.append(0) self.length_required.append(0) for area_i in range(1, self.params['num_tapers']+1): self.area[area_i] = (math.pi / 4) * (self.params['rod_diameter_data'][area_i] ** 2) for i in range(1, self.params['num_tapers']+1): self.a[i] = 1000 * (32.2 * self.params['rod_youngs_modulus_data'][i] * self.area[i] / self.params['rod_weight_data'][i])**(0.5) self.rod_depth[i] = self.rod_depth[i-1] + self.params['rod_length_data'][i] self.pressure[i] = self.pressure[i-1] + self.params['fluid_gradient'] * self.params['rod_length_data'][i] self.buoyant_force[i] = self.pressure[i] * (self.area[i+1] - self.area[i]) self.rod_weight_in_air[i] = self.params['rod_weight_data'][i] * self.params['rod_length_data'][i] self.rod_weight_in_fluid[i] = self.rod_weight_in_air[i] + self.buoyant_force[i] for j in range(1, self.params['num_tapers']+1): for k in range(j+1, self.params['num_tapers']+1): self.weight_data[j] = self.weight_data[j] + self.params['rod_weight_data'][k] * self.params['rod_length_data'][k] for l in range(j, self.params['num_tapers']): self.annular_force_data[j] = self.annular_force_data[j] + self.pressure[l] * (self.area[l] - self.area[l+1]) self.force[j] = (-self.area[self.params['num_tapers']] * self.pressure[self.params['num_tapers']]) + self.weight_data[j] - self.annular_force_data[j] self.alpha[j] = (self.force[j] + self.params['rod_weight_data'][j] * self.params['rod_length_data'][j]) / (self.params['rod_youngs_modulus_data'][j] * 10**6 * self.area[j]) self.stretch[j] = self.stretch[j-1] + self.alpha[j] * self.params['rod_length_data'][j] - (self.params['rod_weight_data'][j] * self.params['rod_length_data'][j]**2) / (2 * self.params['rod_youngs_modulus_data'][j] * 10**6 * self.area[j]) for m in range(1, self.params['num_tapers']+1): self.x_over_a[m] = self.params['rod_length_data'][m] / self.a[m] self.lag_index_array[m] = math.trunc(self.params['rod_length_data'][m] / (self.a[m] * self.dt)) self.factor_array[m] = (self.x_over_a[m] - self.lag_index_array[m] * self.dt) / self.dt self.center_point[m] = self.lag_index_array[m] + 2 self.length_required[m] = 2 * (self.lag_index_array[m] + 1) + 1 self.sum_center_point[1] = self.center_point[1] for n in range(2, self.params['num_tapers']): self.sum_center_point[n] = self.sum_center_point[n-1] + self.center_point[n] - 1 for b in self.buoyant_force: self.buoyant_force_total = self.buoyant_force_total + b for r in self.rod_weight_in_air: self.rod_weight_in_air_total = self.rod_weight_in_air_total + r for r1 in self.rod_weight_in_fluid: self.rod_weight_in_fluid_total = self.rod_weight_in_fluid_total + r1 for r2 in self.params['rod_length_data']: self.rod_depth_total = self.rod_depth_total + r2 for a1 in self.annular_force_data: self.annular_force_data_total = self.annular_force_data_total + a1 for w in self.weight_data: self.weight_data_total = self.weight_data_total + w self.tubing_cross_sectional_area = (math.pi / 4) * (self.params['tubing_od']**2 - self.params['tubing_id']**2) def set_real_dt(self, dt): self.dt = dt for m in range(1, self.params['num_tapers']+1): self.lag_index_array[m] = math.trunc(self.params['rod_length_data'][m] / (self.a[m] * self.dt)) self.factor_array[m] = (self.x_over_a[m] - self.lag_index_array[m] * self.dt) / self.dt def printSetup(self): print("------------ Well Params ---------------------") print("Tubing Head pressure: ", self.params['tubing_head_pressure']) print("Fluid Gradient: ", self.params['fluid_gradient']) print("Stuffing Box Friction: ", self.params['stuffing_box_friction']) print("Number of Tapers: ", self.params['num_tapers']) print("Damping Factor: " , self.params['c'] ) print("Rod Length Data: " , self.params['rod_length_data']) print("Rod Diameter Data: ", self.params['rod_diameter_data']) print("Rod Young's Modulus Data: ", self.params['rod_youngs_modulus_data']) print("Rod Weight Data: ", self.params['rod_weight_data']) print("dt: ", self.dt) print("tubing_id:", self.params['tubing_id']) print("tubing_od:", self.params['tubing_od']) print("----------------------------------------------") def printTaper(self): print("------------ Taper Params --------------------") print("area:", self.area) print("Speed of Sound in Rod(a):", self.a) print("Rod Depth:", self.rod_depth) print("pressure:", self.pressure) print("Buoyant force:", self.buoyant_force) print("Rod Weight in Air:", self.rod_weight_in_air) print("Weight Data:", self.weight_data) print("Annulus pressure Data:", self.annular_force_data) print("force Data:", self.force) print("alpha Data:", self.alpha) print("stretch Data:", self.stretch) print("X over A:", self.x_over_a) print("Factor Array:", self.factor_array) print("Lag Index Array:", self.lag_index_array) print("center_point:", self.center_point) print("sum_center_point:", self.sum_center_point) print("Length Required:", self.length_required) print("TubingCSA:", self.tubing_cross_sectional_area) print("------------ Taper Totals --------------------") print("buoyant_force_total:", self.buoyant_force_total) print("rod_weight_in_air_total:", self.rod_weight_in_air_total) print("rod_weight_in_fluid_total:", self.rod_weight_in_fluid_total) print("rod_depth_total:", self.rod_depth_total) print("annular_force_data_total:", self.annular_force_data_total) print("weight_data_total:", self.weight_data_total) print("----------------------------------------------") def position_function(p,tap): #CHECKED LOGIC ON 1/21/2015 global top_load_array global top_pos_array global load_before global load_after global load_before_3 global load_after_3 load_before = 0 load_after = 0 load_before_3 = 0 load_after_3 = 0 dt = tap.dt a1= tap.a[p] cd = tap.params['c'][p] factor = tap.factor_array[p] rod_length = tap.params['rod_length_data'][p] lag_index = tap.lag_index_array[p] Y1 = tap.params['rod_youngs_modulus_data'][p] Y1 = Y1 * 10**6 A1 = tap.area[p] center_of_array = tap.center_point[p] #print("p:", p, "a1:", a1, "| cd:", cd, "| factor:", factor, "| rod_length:", rod_length, "| lag_index:", lag_index, "| Y1:", Y1, "| A1:", A1, "| center_of_array:", center_of_array) i_before = center_of_array - lag_index i_after = center_of_array + lag_index #print("------") pump_pos = 0 pump_pos = math.exp((cd * rod_length) / (2 * a1)) * (top_pos_array[p][i_after] + factor * (top_pos_array[p][i_after + 1] - top_pos_array[p][i_after])) pump_pos = pump_pos + math.exp(-(cd * rod_length) / (2 * a1)) * (top_pos_array[p][i_before] + factor * (top_pos_array[p][i_before - 1] - top_pos_array[p][i_before])) pump_pos = 0.5 * pump_pos inside_integral = 0 q = 1 while(q < 2* lag_index - 1): inside_integral = inside_integral + dt * (1 / (Y1 * A1)) * (math.exp(-(cd * (lag_index - q ) * dt) / 2 ) * top_load_array[p][i_before +q]) q+=1 inside_integral = inside_integral + (1/2) * dt * (1 / (Y1 * A1)) * (math.exp(-(cd *lag_index * dt)/ 2) * top_load_array[p][i_before] + math.exp(- (cd* (-lag_index ) * dt) / 2) * top_load_array[p][i_after]) load_before = math.exp(-(cd * lag_index * dt) / 2) * top_load_array[p][i_before] + factor * (math.exp(-(cd * (lag_index + 1) * dt) / 2) * top_load_array[p][i_before - 1] - math.exp(-(cd* (lag_index) * dt ) / 2 ) * top_load_array[p][i_before]) load_after = math.exp(-(cd * (-lag_index) * dt) / 2) * top_load_array[p][i_after] + factor * (math.exp(-(cd * (-lag_index - 1) * dt) / 2) * top_load_array[p][i_after + 1] - math.exp(-(cd* (-lag_index) * dt ) / 2 ) * top_load_array[p][i_after]) inside_integral = inside_integral + 0.5 * factor * dt * (1 / (Y1 * A1)) * (load_before + math.exp(-(cd * (lag_index) * dt ) / 2) * top_load_array[p][i_before]) inside_integral = inside_integral + 0.5 * factor * dt * (1 / (Y1 * A1)) * (load_after + math.exp(-(cd * (-lag_index) * dt ) / 2) * top_load_array[p][i_after]) inside_integral = (1/2) * a1 * inside_integral pump_pos = pump_pos + inside_integral inside_integral = 0 r = 1 while(r < 2 * lag_index - 1): inside_integral = inside_integral + dt * (math.exp(-(cd *(lag_index - r) * dt)/2) * (top_pos_array[p][i_before + r])) r = r+1 inside_integral = inside_integral + (1/2) * dt * (math.exp(-(cd * lag_index * dt) / 2) * top_pos_array[p][i_before] + math.exp(-(cd * (-lag_index) * dt)/2) * top_pos_array[p][i_after]) load_before_3 = math.exp(-(cd * (lag_index) *dt) / 2) * top_pos_array[p][i_before] + factor * (math.exp(-(cd * (lag_index + 1) *dt) / 2) * top_pos_array[p][i_before - 1] - math.exp(-(cd * (lag_index) * dt)/2) * top_pos_array[p][i_before]) load_after_3 = math.exp(-(cd * (-lag_index) *dt) / 2) * top_pos_array[p][i_after] + factor * (math.exp(-(cd * (-lag_index - 1) *dt) / 2) * top_pos_array[p][i_after + 1] - math.exp(-(cd * (-lag_index) * dt)/2) * top_pos_array[p][i_after]) inside_integral = inside_integral + (1/2) * factor * dt * (load_before_3 + math.exp(-(cd * (lag_index) * dt)/2) * top_pos_array[p][i_before]) inside_integral = inside_integral + (1/2) * factor * dt * (load_after_3 + math.exp(-(cd * (-lag_index) * dt)/2) * top_pos_array[p][i_after]) inside_integral = -((cd * rod_length) / 4 ) * ((1/2) * (cd / (2*a1))) * inside_integral pump_pos = pump_pos + inside_integral return(pump_pos) def load_function(s, tap): global top_load_array global top_pos_array global load_before global load_after global load_before_3 global load_after_3 dt = tap.dt a1= tap.a[s] cd = tap.params['c'][s] factor = tap.factor_array[s] rod_length = tap.params['rod_length_data'][s] lag_index = tap.lag_index_array[s] Y1 = tap.params['rod_youngs_modulus_data'][s] Y1 = Y1 * 10**6 A1 = tap.area[s] center_of_array = tap.center_point[s] i_before = center_of_array - lag_index i_after = center_of_array + lag_index pump_load = 0 pump_load = (1/2) * (a1 / (Y1*A1)) * (1 / a1) * (load_before + load_after) temp_result = Y1 * A1 * pump_load pump_load = pump_load -((cd * rod_length) / 4 ) * ((1/2) * (cd / (2*a1))) * (1/a1) * (load_before_3 + load_after_3) first_part = 0 point_after = (top_pos_array[s][i_after+1] - top_pos_array[s][i_after -1]) / (2 * dt) point_before = (top_pos_array[s][i_before+1] - top_pos_array[s][i_before -1]) / (2 * dt) first_part = (math.exp((cd*rod_length)/(2*a1)) * point_after - math.exp(-(cd*rod_length)/(2*a1)) * point_before) / (2 *a1) first_part = first_part +(cd * math.exp((cd*rod_length)/(2*a1))*top_pos_array[s][i_after] - cd * math.exp((-cd*rod_length)/(2*a1))*top_pos_array[s][i_before]) / (4 * a1) pump_load = Y1 * A1 * (first_part + pump_load) return(pump_load) def calc(pos, load, t): global pr_p global pr_l global pmp_p global pmp_l global top_pos_array global top_load_array global position_load_counter global surface_pos_array global surface_load_array USE_SHIFT = False surface_pos_array.append(pos) surface_load_array.append(load) load_mult = 1 tapers_allowed = 1 pr_l[0] = pr_l[1] pmp_l[0] = pmp_l[1] pr_p[0] = pr_p[1] pmp_p[0] = pmp_p[1] pr_p[1] = pos pr_l[1] = load for ii in range(1, t.length_required[1]+1): top_pos_array[1][ii-1] = top_pos_array[1][ii] top_load_array[1][ii-1] = top_load_array[1][ii] top_pos_array[1][t.length_required[1]] = -(pr_p[1] / 12) if(pr_p[1] > pr_p[0]): top_load_array[1][t.length_required[1]] = load_mult * (pr_l[1] - t.rod_weight_in_fluid_total) - t.params['stuffing_box_friction'] elif (pr_p[1] < pr_p[0]): top_load_array[1][t.length_required[1]] = load_mult * (pr_l[1] - t.rod_weight_in_fluid_total) + t.params['stuffing_box_friction'] else: top_load_array[1][t.length_required[1]] = load_mult * (pr_l[1] - t.rod_weight_in_fluid_total) j=1 while j <= tapers_allowed: count[j] = count[j] + 1 if (count[j] >= t.length_required[j]): if((j+1) <= t.params['num_tapers']): for ii in range(2, t.length_required[j+1]+1): top_pos_array[j+1][ii-1] = top_pos_array[j+1][ii] top_load_array[j+1][ii-1] = top_load_array[j+1][ii] top_pos_array[j+1][t.length_required[j+1]] = position_function(j,t) top_load_array[j+1][t.length_required[j+1]] = load_function(j,t) else: if USE_SHIFT: pmp_p[1] = -12 * (position_function(j,t) + t.stretch[t.params['num_tapers']]) else: pmp_p[1] = -12 * position_function(j,t) pmp_l[1] = load_function(j,t) + t.force[t.params['num_tapers']] pump_pos_array.append(pmp_p[1]) pump_load_array.append(pmp_l[1]) position_load_counter += 1 count[j] = count[j] - 1 tapers_allowed = tapers_allowed + 1 if(tapers_allowed > t.params['num_tapers']): tapers_allowed = t.params['num_tapers'] j = j+1 return [[pr_p[1], pr_l[1]], [pmp_p[1], pmp_l[1]]] tapers = Taper(taper_setup) if __name__ == '__main__': file_path = fd.askopenfilename() store_points = False import csv with open(file_path, 'r') as csvfile: data_reader = csv.reader(csvfile) for row in data_reader: if (row[0]=="s_pos"): store_points = False if store_points: surface.append([float(row[0]),float(row[1])]) if (row[0]=="d_pos"): store_points = True def loop(iterations): stroke_start_time = time.time() point_times = [0,0] for point in range(0,len(surface)): #print ("calc for", surface[point][0], surface[point][1] ) point_time_delta = point_times[1] - point_times[0] #print("pt took:", point_time_delta) point_times[0] = point_times[1] point_times[1] = time.time() current = calc(surface[point][0], surface[point][1], tapers) if ((point_time_delta > 0) and (point_time_delta < 2)): tapers.set_real_dt(point_time_delta) time.sleep(tapers.params['dt']) stroke_end_time = time.time() stroke_time_delta = stroke_end_time - stroke_start_time pump_pos_array.append(pump_pos_array[0]) pump_load_array.append(pump_load_array[0]) surface_pos_array.append(surface_pos_array[0]) surface_load_array.append(surface_load_array[0]) last_card = Card(1,surface_pos_array, surface_load_array, pump_pos_array, pump_load_array, tapers, unit_config) last_card.set_strokes_per_minute((1 / stroke_time_delta) * 60) last_card.calc_fillage(25,25) last_card.calc_hp(100) print("polished_rod_horsepower:", last_card.polished_rod_horsepower) print("pump_horsepower:", last_card.pump_horsepower) print("strokes_per_minute:", last_card.strokes_per_minute) iterations -= 1 if (iterations > 0): loop(iterations) if (iterations == 0): last_card.print_corners() last_card.write_csv() last_card.plot() loop(1) #pprint(vars(last_card))