Files
POC-Python/algorithm.py
Patrick McDonagh 0adef36b26 Adds files
2016-05-27 17:38:12 -05:00

727 lines
34 KiB
Python

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))