import math import matplotlib matplotlib.use("TkAgg") import matplotlib.pyplot as plt from collections import deque import datetime import time import csv import json import setupTapers print("Be sure that this folder includes setupTapers.py and that you've put a surface card in surface_card.json.") ready = str(raw_input("Are you ready to start? (y/n): ")).lower() if not (ready == 'y'): raise Exception("Not ready to start!") well_setup = setupTapers.get_setup() surface = [] with open('surface_card.json') as sc: surface = json.loads(sc.read()) load_before = 0.0 load_after = 0.0 load_before_3 = 0.0 load_after_3 = 0.0 blank_array = [0.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, 0.0] pr_p = [0.0, 0.0] pmp_p = [0.0, 0.0] pmp_l = [0.0, 0.0] count = [0] * (well_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.well_setup = 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") self.riemann_slices = 100 def set_strokes_per_minute(self, strokes_per_minute): self.strokes_per_minute = strokes_per_minute 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 distance_to_line(self, x0, y0): """ Finds the perpendicular distance from a point to the line between max and min points""" x1 = self.downhole_min_position[0] x2 = self.downhole_max_position[0] y1 = self.downhole_min_load[1] y2 = self.downhole_max_load[1] d = abs((y2-y1)*x0 - (x2-x1)*y0 + x2*y1 - y2*x1) / math.sqrt(math.pow(y2-y1, 2) + math.pow(x2-x1, 2)) return d def find_corners(self, sl=100): self.find_limits() def lineresolve(x1, x2, y1, y2, targ): m = ((y2 - y1) / (x2 - x1)) b = y1 - m * x1 return(m * targ + b) slices = int(sl) num_points = len(self.downhole_pos) dh_pos = self.downhole_pos dh_lod = self.downhole_load delta_p = (self.downhole_max_position[0] - self.downhole_min_position[0]) / float(slices) position_targets = [] for i in range(0, slices): position_targets.append(self.downhole_min_position[0] + i * delta_p) self.top_slices = [] last_pos_index = 0 for i in range(0, slices): targ = position_targets[i] for j in range(last_pos_index, num_points - 1): if (dh_pos[j] <= targ and dh_pos[j+1] > targ): found_i = j next_i = j+1 fake_load = lineresolve(dh_pos[found_i], dh_pos[next_i], dh_lod[found_i], dh_lod[next_i], targ) self.top_slices.append((targ, fake_load)) last_pos_index = found_i break self.bot_slices = [] last_pos_index = 0 for i in range(0, slices): targ = position_targets[i] for j in range(0, num_points - 1): if (dh_pos[j] > targ and dh_pos[j+1] <= targ): found_i = j next_i = j+1 fake_load = lineresolve(dh_pos[found_i], dh_pos[next_i], dh_lod[found_i], dh_lod[next_i], targ) self.bot_slices.append((targ, fake_load)) break top_d = [] bot_d = [] # Here's where we get the distance from each point to the line for i in range(1, slices): try: top_d.append(self.distance_to_line(self.top_slices[i][0], self.top_slices[i][1])) except Exception as e: print e print("Error - top_d, index: {0}".format(i)) try: bot_d.append(self.distance_to_line(self.bot_slices[i][0], self.bot_slices[i][1])) except Exception: print("Error - bot_d, index: {0}".format(i)) top_cor_i = top_d.index(max(top_d)) bot_cor_i = bot_d.index(max(bot_d)) self.corner_top_left = self.top_slices[top_cor_i] self.corner_fillage = self.bot_slices[bot_cor_i] self.corner_top_right = self.downhole_max_position self.corner_bottom_left = self.downhole_min_position def print_corners(self): if self.corner_fillage: self.find_corners(100) print("Bottom Left:", self.corner_bottom_left) print("Top Left:", self.corner_top_left) print("Top Right:", self.corner_top_right) print("Fill Point:", self.corner_fillage) def calc_fillage(self): if not ('self.corner_fillage' in locals()): self.find_corners(100) self.downhole_gross_stroke = self.downhole_max_position[0] - self.downhole_min_position[0] self.tubing_movement = 12 * (tapers.rod_depth_total - self.well_setup['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_and_analyze(self): self.polished_rod_horsepower = 0 self.pump_horsepower = 0 if not self.surface_stroke_length: self.find_limits() if not self.downhole_net_stroke: self.calc_fillage() dx_surface = self.surface_stroke_length / (self.riemann_slices + 1) dx_downhole = self.downhole_net_stroke / (self.riemann_slices + 1) for i in range(1, self.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) top_slices.push(slice_downhole[1]) bottom_slices.push(slice_downhole[2]) self.polished_rod_horsepower += (dx_surface / 12) * slice_surface[0] * self.strokes_per_minute / 33000 self.pump_horsepower += (dx_downhole / 12) * slice_downhole[0] * self.strokes_per_minute / 33000 # CARD ANALYSIS / $$$ MAKER def prepare_points(analysis_slices): if (analysis_slices > (self.riemann_slices * 0.50)): return 911 top_left_average = 0 top_right_average = 0 top_mid_average = 0 bot_left_average = 0 bot_right_average = 0 bot_mid_average = 0 pct_slices = (analysis_slices / self.riemann_slices) * 100 first_x_pct_end = math.floor(self.riemann_slices*pct_slices) x_pct_ct = math.floor(self.riemann_slices*pct_slices) last_x_pct_start = self.riemann_slices - math.floor(self.riemann_slices * pct_slices) num_btw_slices = self.riemann_slices - (3 * x_pct_ct) # number of slices between the measurements used for top/bottom average and the mid point average # Find average of leftmost Points for la_i in range(0, first_x_pct_end): top_left_average += self.top_slices[la_i] * (1/x_pct_ct) bot_left_average += self.bottom_slices[la_i] * (1/x_pct_ct) # Find average of rightmost points for ra_i in range(last_x_pct_start, self.riemann_slices): top_right_average += self.top_slices[ra_i] * (1/x_pct_ct) bot_right_average += self.bottom_slices[ra_i] * (1/x_pct_ct) # Find average of middle points for ma_i in range(math.floor((self.riemann_slices/2) - self.riemann_slices * (pct_corners / 2)), math.floor((self.riemann_slices / 2) + self.riemann_slices * (pct_corners / 2))): top_mid_average += self.top_slices[ma_i] * (1/x_pct_ct) bot_mid_average += self.bottom_slices[ma_i] * (1/x_pct_ct) 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.well_setup['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_and_analyze(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.well_setup['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') 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.0 self.rod_weight_in_fluid_total = 0.0 self.rod_weight_in_air_total = 0.0 self.weight_data_total = 0.0 self.rod_depth_total = 0.0 self.annular_force_data_total = 0.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.0) self.pressure.append(0.0) self.buoyant_force.append(0.0) self.rod_weight_in_air.append(0.0) self.rod_depth.append(0.0) self.rod_weight_in_fluid.append(0.0) self.weight_data.append(0.0) self.annular_force_data.append(0.0) self.a.append(0.0) self.stretch.append(0.0) self.force.append(0.0) self.alpha.append(0.0) self.x_over_a.append(0.0) self.factor_array.append(0.0) self.lag_index_array.append(0.0) self.center_point.append(0.0) self.sum_center_point.append(0.0) self.length_required.append(0.0) for area_i in range(1, self.params['num_tapers']+1): self.area[area_i] = (math.pi / 4.0) * (self.params['rod_diameter_data'][area_i] ** 2.0) for i in range(1, self.params['num_tapers']+1): self.a[i] = 1000.0 * (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.0 load_after = 0.0 load_before_3 = 0.0 load_after_3 = 0.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.0 pump_pos = math.exp((cd * rod_length) / (2.0 * 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.0 * 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.0 q = 1 while(q < 2 * lag_index - 1): inside_integral = inside_integral + dt * (1.0 / (Y1 * A1)) * (math.exp(-(cd * (lag_index - q) * dt) / 2.0) * top_load_array[p][i_before + q]) q += 1 inside_integral = inside_integral + (0.5) * dt * (1.0 / (Y1 * A1)) * (math.exp(-(cd * lag_index * dt) / 2.0) * top_load_array[p][i_before] + math.exp(- (cd * (-lag_index) * dt) / 2.0) * top_load_array[p][i_after]) load_before = math.exp(-(cd * lag_index * dt) / 2.0) * top_load_array[p][i_before] + factor * (math.exp(-(cd * (lag_index + 1.0) * dt) / 2.0) * top_load_array[p][i_before - 1] - math.exp(-(cd * (lag_index) * dt) / 2.0) * top_load_array[p][i_before]) load_after = math.exp(-(cd * (-lag_index) * dt) / 2.0) * top_load_array[p][i_after] + factor * (math.exp(-(cd * (-lag_index - 1.0) * dt) / 2.0) * top_load_array[p][i_after + 1] - math.exp(-(cd * (-lag_index) * dt) / 2.0) * top_load_array[p][i_after]) inside_integral = inside_integral + 0.5 * factor * dt * (1.0 / (Y1 * A1)) * (load_before + math.exp(-(cd * (lag_index) * dt) / 2.0) * top_load_array[p][i_before]) inside_integral = inside_integral + 0.5 * factor * dt * (1.0 / (Y1 * A1)) * (load_after + math.exp(-(cd * (-lag_index) * dt) / 2.0) * top_load_array[p][i_after]) inside_integral = (0.5) * a1 * inside_integral pump_pos = pump_pos + inside_integral inside_integral = 0.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 + (0.5) * dt * (math.exp(-(cd * lag_index * dt) / 2.0) * top_pos_array[p][i_before] + math.exp(-(cd * (-lag_index) * dt)/2.0) * top_pos_array[p][i_after]) load_before_3 = math.exp(-(cd * (lag_index) * dt) / 2.0) * top_pos_array[p][i_before] + factor * (math.exp(-(cd * (lag_index + 1) * dt) / 2.0) * top_pos_array[p][i_before - 1] - math.exp(-(cd * (lag_index) * dt) / 2.0) * top_pos_array[p][i_before]) load_after_3 = math.exp(-(cd * (-lag_index) * dt) / 2.0) * top_pos_array[p][i_after] + factor * (math.exp(-(cd * (-lag_index - 1.0) * dt) / 2.0) * top_pos_array[p][i_after + 1] - math.exp(-(cd * (-lag_index) * dt) / 2.0) * top_pos_array[p][i_after]) inside_integral = inside_integral + (0.5) * factor * dt * (load_before_3 + math.exp(-(cd * (lag_index) * dt) / 2.0) * top_pos_array[p][i_before]) inside_integral = inside_integral + (0.5) * factor * dt * (load_after_3 + math.exp(-(cd * (-lag_index) * dt) / 2.0) * top_pos_array[p][i_after]) inside_integral = -((cd * rod_length) / 4.0) * ((0.5) * (cd / (2.0 * 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] 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.0**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.0 pump_load = (0.5) * (a1 / (Y1*A1)) * (1.0 / a1) * (load_before + load_after) pump_load = pump_load - ((cd * rod_length) / 4.0) * ((0.5) * (cd / (2.0 * a1))) * (1.0 / a1) * (load_before_3 + load_after_3) first_part = 0.0 point_after = (top_pos_array[s][i_after+1] - top_pos_array[s][i_after - 1]) / (2.0 * dt) point_before = (top_pos_array[s][i_before+1] - top_pos_array[s][i_before - 1]) / (2.0 * dt) first_part = (math.exp((cd*rod_length)/(2.0*a1)) * point_after - math.exp(-(cd*rod_length)/(2*a1)) * point_before) / (2.0 * a1) first_part = first_part + (cd * math.exp((cd*rod_length)/(2.0*a1))*top_pos_array[s][i_after] - cd * math.exp((-cd*rod_length)/(2*a1))*top_pos_array[s][i_before]) / (4.0 * 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.0) 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(well_setup) if __name__ == '__main__': 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) print(current) 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, well_setup) last_card.set_strokes_per_minute((1 / stroke_time_delta) * 60.0) last_card.calc_fillage() last_card.print_corners() # last_card.calc_hp_and_analyze(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(2) # pprint(vars(last_card))