Contents
- Constants (SHOULD THESE BE ELSEWHERE?)
- Instrumentation
- Current Values
- User Inputs
- Production Parameters
- Rod String Inputs
- Calculated Taper Parameters
- Constructor Function
- Well Setup Function
- Position Function (From Algorithm)
- Load Function (From Algorithm)
- Calc Function
- Check if the stroke has ended
- Evaluation function for simulated position and load points
- End of Stroke actions
- Draws Cards
- Determines Rod Weight Per Foot given Young's Modulus and Diameter
classdef Well < handle + properties(Constant) +
Constants (SHOULD THESE BE ELSEWHERE?)
YM_STEEL=30.5; + YM_FIBERGLASS=7.2; + + DIRECTION_UNKNOWN=0; + DIRECTION_UP=1; + DIRECTION_DOWN=2; +
end
+
+ properties
+Instrumentation
device; + mux; + positionSensor; + loadSensor; +
Current Values
currentSurfacePosition,currentSurfaceLoad; + currentDownholePosition,currentDownholeLoad; + currentCard=Card(0); + lastCard=Card(0); + direction=0; + lastDirection=0; +
User Inputs
wellName; + dt = 1.0 / 30.0; + tubingHeadPressure = 40.0; + fluidGradient = 0.45; + stuffingBoxFriction = 100.0; + numTapers = 3; + tubingAnchorDepth = 0.0; + pumpDiameter = 2.5; + pumpArea; + tubingInnerDiameter = 1.995; + tubingOuterDiameter = 2.375; + tubingCrossSectionalArea; + structuralRating = 320000; + frictionEstimate; +
Production Parameters
kFactor=1.0; + waterBblRatio=0.90; + oilBblRatio=0.10; + gasMcfRatio=5.0; +
Rod String Inputs
c; + rodLength; + rodDiameter; + rodYM; + rodWeightPerFoot; +
Calculated Taper Parameters
a; + area; + pressure; + stretch; + force; + alpha; + xOverA; + factorArray; + lagIndex; + lengthRequired; + centerPoint; + count; + + buoyantForceTotal = 0.0; + rodDepthTotal = 0.0; + rodWeightAirTotal = 0.0; + rodWeightFluidTotal = 0.0; +
end + + properties(Access=private) + topPosArray; + topLoadArray; + + % cross-function variables for position and load + % POSITION MUST ALWAYS BE CALLED FIRST. + loadBefore = 0.0; + loadAfter = 0.0; + loadBefore3 = 0.0; + loadAfter3 = 0.0; + tempLoadValue = 0.0; + end + + methods +
Constructor Function
function me = Well(name, sim) + me.wellName = name; + + me.c = [0.8 0.8 0.8]; + me.rodLength = [1475.0 1525.0 2000.0]; + me.rodDiameter = [1.0 0.875 0.750]; + me.rodYM = [me.YM_STEEL me.YM_STEEL me.YM_STEEL]; + me.rodWeightPerFoot = [Well.lookupRodWeightPerFoot(me.rodYM(1), me.rodDiameter(1)) ... + Well.lookupRodWeightPerFoot(me.rodYM(2), me.rodDiameter(2)) ... + Well.lookupRodWeightPerFoot(me.rodYM(3), me.rodDiameter(3))]; + + % Initialize + me.a = zeros(1, me.numTapers, 'double'); + me.area = zeros(1, me.numTapers+1, 'double'); + me.pressure = zeros(1, me.numTapers, 'double'); + me.stretch = zeros(1, me.numTapers, 'double'); + me.force = zeros(1, me.numTapers, 'double'); + me.alpha = zeros(1, me.numTapers, 'double'); + me.xOverA = zeros(1, me.numTapers, 'double'); + me.factorArray = zeros(1, me.numTapers, 'double'); + me.lagIndex = zeros(1, me.numTapers, 'uint8'); + me.lengthRequired = zeros(1, me.numTapers, 'uint8'); + me.centerPoint = zeros(1, me.numTapers, 'uint8'); + me.count = zeros(1, me.numTapers, 'uint8'); + me.topPosArray = zeros(10, 250,'double'); + me.topLoadArray = zeros(10, 250, 'double'); + + if logical(sim) + me.device = 0; + me.mux = MuxSetupSim(me.device); + me.positionSensor = AnalogInputSim(me.mux, 1, 0.0, 65535.0, 0.0, 65535.0); + me.loadSensor = AnalogInputSim(me.mux, 2, 0.0, 65535.0, 0.0, 65535.0); + else + me.device = raspi('10.0.0.106', 'pi', 'HenryPump@1903'); + me.mux = MuxSetupSim(me.device); + me.positionSensor = AnalogInput(mux, 1, 0.0, 65535.0, 0.0, 140.0); + me.loadSensor = AnalogInput(mux, 2, 0.0, 65535.0, 0.0, 50000.0); + end + + me.wellSetup(); + end +
Not enough input arguments. + +Error in Well (line 96) + me.wellName = name; +
Well Setup Function
function wellSetup(me) + me.tubingCrossSectionalArea = (pi / 4) * ... + (power(me.tubingOuterDiameter, 2) - ... + power(me.tubingInnerDiameter,2)); + me.pumpArea = power(me.pumpDiameter, 2) * pi; + + for i = 1:me.numTapers + me.area(i) = pi / 4 * power(me.rodDiameter(i), 2); + end + + for i = 1:me.numTapers + me.a(i) = 1000.0 * sqrt(32.2 * me.rodYM(i) * me.area(i) / me.rodWeightPerFoot(i)); + me.rodDepthTotal = me.rodDepthTotal + me.rodLength(i); + + if i > 1 + me.pressure(i) = me.pressure(i - 1) + me.fluidGradient * me.rodLength(i); + else + me.pressure(i) = me.tubingHeadPressure + me.fluidGradient * me.rodLength(i); + end + + me.buoyantForceTotal = me.buoyantForceTotal + me.pressure(i) * (me.area(i+1) - me.area(i)); + me.rodWeightAirTotal = me.rodWeightAirTotal + me.rodWeightPerFoot(i) * me.rodLength(i); + me.rodWeightFluidTotal = me.rodWeightAirTotal + me.buoyantForceTotal; + end + + for j = 1:me.numTapers + weightData = 0.0; + annularForceData = 0.0; + + for i = j+1:me.numTapers + weightData = weightData + me.rodWeightPerFoot(i) * me.rodLength(i); + end + + for i = j:me.numTapers-1 + annularForceData = annularForceData + me.pressure(i) * (me.area(i) - me.area(i+1)); + end + + me.force(j) = (-me.area(me.numTapers) * me.pressure(me.numTapers)) + weightData - annularForceData; + me.alpha(j) = (me.force(j) + me.rodWeightPerFoot(j) * me.rodLength(j)) / (me.rodYM(j) * power(10, 6) * me.area(j)); + + if j > 1 + me.stretch(j) = me.stretch(j - 1) + me.alpha(j) * me.rodLength(j) - ... + (me.rodWeightPerFoot(j) * power(me.rodLength(j), 2.0)) / ... + (2.0 * me.rodYM(j) * power(10, 6) * me.area(j)); + else + me.stretch(j) = 0.0 + me.alpha(j) * me.rodLength(j) - ... + (me.rodWeightPerFoot(j) * power(me.rodLength(j), 2.0)) / ... + (2.0 * me.rodYM(j) * power(10, 6) * me.area(j)); + end + end + + for i = 1:me.numTapers + me.xOverA(i) = me.rodLength(i) / me.a(i); + me.lagIndex(i) = floor(me.rodLength(i) / (me.a(i) * me.dt)); + me.factorArray(i) = (me.xOverA(i) - double(me.lagIndex(i)) * me.dt) / me.dt; + me.centerPoint(i) = me.lagIndex(i) + 2; + me.lengthRequired(i) = 2 * (me.lagIndex(i) + 1) + 1; + end + + me.frictionEstimate = me.rodDepthTotal * 0.10; + end +
Position Function (From Algorithm)
function pumpPosition = position(me, i) + % Position Function + ai = me.a(i); + ci = me.c(i); + factori = me.factorArray(i); + lengthi = me.rodLength(i); + lagi = double(me.lagIndex(i)); + yi = me.rodYM(i) * power(10, 6); + areai = me.area(i); + centeri = me.centerPoint(i); + iBefore = centeri - me.lagIndex(i); + iAfter = centeri + me.lagIndex(i); + + me.loadBefore = 0.0; + me.loadAfter = 0.0; + me.loadBefore3 = 0.0; + me.loadAfter3 = 0.0; + + pumpPosition = exp(ci * lengthi / (2.0 * ai)) * (me.topPosArray(i, iAfter) + ... + factori * (me.topPosArray(i, iAfter+1) - me.topPosArray(i, iAfter))); + + pumpPosition = pumpPosition + exp(-ci * lengthi / (2.0 * ai)) * ... + (me.topPosArray(i,iBefore) + factori * (me.topPosArray(i,iBefore-1) - ... + me.topPosArray(i, iBefore))); + + pumpPosition = 0.5 * pumpPosition; + insideIntegral = 0.0; + + for jj = 1:2*lagi-1 + insideIntegral = insideIntegral + me.dt / (yi * areai) * ... + (exp((-ci * (lagi - jj) *me.dt) / 2.0) * ... + me.topLoadArray(i, iBefore + jj)); + end + + insideIntegral = insideIntegral + 0.5 * me.dt /(yi * areai) * ... + (exp((-ci * lagi * me.dt) / 2.0) * me.topLoadArray(i, iBefore) + ... + exp((-ci * -lagi * me.dt) / 2.0) * me.topLoadArray(i, iAfter)); + + me.loadBefore = exp((-ci * lagi * me.dt) / 2.0) * me.topLoadArray(i, iBefore) + ... + factori * (exp((-ci * (lagi + 1) * me.dt) / 2.0) * me.topLoadArray(i, iBefore-1) - ... + exp((-ci * lagi * me.dt) / 2.0) * me.topLoadArray(i, iBefore)); + + me.loadAfter = exp((-ci * -lagi * me.dt) / 2.0) * me.topLoadArray(i, iAfter) + ... + factori * (exp((-ci * (-lagi - 1) * me.dt) / 2.0) * me.topLoadArray(i, iAfter+1) - ... + exp((-ci * -lagi * me.dt) / 2.0) * me.topLoadArray(i,iAfter)); + + insideIntegral = insideIntegral + 0.5 * factori * me.dt / (yi * areai) * ... + (me.loadBefore + exp((-ci * lagi * me.dt) / 2.0) * me.topLoadArray(i, iBefore)); + + insideIntegral = insideIntegral + 0.5 * factori * me.dt / (yi * areai) * ... + (me.loadAfter + exp((-ci * -lagi * me.dt) / 2.0) * me.topLoadArray(i, iAfter)); + + insideIntegral = 0.5 * ai * insideIntegral; + pumpPosition = pumpPosition + insideIntegral; + + insideIntegral = 0.0; + + for jj = 1:2*lagi-1 + insideIntegral = insideIntegral + me.dt * (exp((-ci * (lagi - jj) * me.dt) / 2.0) * ... + me.topPosArray(i, iBefore+jj)); + end + + insideIntegral = insideIntegral + 0.5 * me.dt * (exp((-ci * lagi * me.dt) / 2.0) * ... + me.topPosArray(i, iBefore) + exp((-ci * -lagi * me.dt) / 2.0) * ... + me.topPosArray(i, iAfter)); + + me.loadBefore3 = exp((-ci * lagi * me.dt) / 2.0) * me.topPosArray(i, iBefore) + ... + factori * (exp((-ci * (lagi+1) * me.dt) / 2.0) * me.topPosArray(i, iBefore-1) - ... + exp((-ci * lagi * me.dt) / 2.0) * me.topPosArray(i, iBefore)); + + me.loadAfter3 = exp((-ci * -lagi * me.dt) / 2.0) * me.topPosArray(i, iAfter) + ... + factori * (exp((-ci * (-lagi-1) * me.dt) / 2.0) * me.topPosArray(i, iAfter+1) - ... + exp((-ci * -lagi * me.dt) / 2.0) * me.topPosArray(i, iAfter)); + + insideIntegral = insideIntegral + 0.5 * factori * me.dt * (me.loadBefore3 + ... + exp((-ci * lagi * me.dt) / 2.0) * me.topPosArray(i, iBefore)); + + insideIntegral = insideIntegral + 0.5 * factori * me.dt * (me.loadAfter3 + ... + exp((-ci * -lagi * me.dt) / 2.0) * me.topPosArray(i, iAfter)); + + insideIntegral = -(ci * lengthi / 4) * 0.5 * (ci / (2.0 * ai)) * insideIntegral; + me.tempLoadValue = insideIntegral / lengthi; + pumpPosition = pumpPosition + insideIntegral; + end +
Load Function (From Algorithm)
function pumpLoad = load(me, i) + % Load Function + ai = me.a(i); + ci = me.c(i); + lengthi = me.rodLength(i); + lagi = me.lagIndex(i); + yi = me.rodYM(i) * power(10, 6); + areai = me.area(i); + centeri = me.centerPoint(i); + iBefore = centeri - lagi; + iAfter = centeri + lagi; + + pumpLoad = 0.5 * (ai / (yi * areai)) * (1 / ai) * (me.loadBefore + me.loadAfter); + tempResult = yi * areai * pumpLoad; + pumpLoad = pumpLoad - (ci * lengthi / 4) * (0.5 * (ci / (2.0 * ai))) * (1 / ai) * ... + (me.loadBefore3 + me.loadAfter3); + + ptAfter = (me.topPosArray(i, iAfter+1) - me.topPosArray(i, iAfter-1)) / (2.0 * me.dt); + ptBefore = (me.topPosArray(i, iBefore+1) - me.topPosArray(i, iBefore-1)) / (2.0 * me.dt); + firstPart = (exp((ci * lengthi) / (2.0 * ai)) * ptAfter - exp((-ci * lengthi) / ... + (2.0 * ai)) * ptBefore) / (2.0 * ai); + + firstPart = firstPart + (ci * exp((ci * lengthi) / (2.0 * ai)) * ... + me.topPosArray(i, iAfter) - ci * exp((-ci * lengthi) / (2.0 * ai)) * ... + me.topPosArray(i, iBefore)) / (4 * ai); + + pumpLoad = yi * areai * (firstPart + pumpLoad); % + tempLoadValue ? + end +
Calc Function
Calls position and load functions
function [pumpPosition, pumpLoad, status] = calc(me, polishedRodPosition, lastPolishedRodPosition,... + polishedRodLoad) + me.currentSurfacePosition = polishedRodPosition; + me.currentSurfaceLoad = polishedRodLoad; + pumpPosition = -1.0; + pumpLoad = -1.0; + status = -1; + loadMult = 1.0; + tapersAllowed = 1; + + for ii = 2:me.lengthRequired(1) + me.topPosArray(1, ii-1) = me.topPosArray(1, ii); + me.topLoadArray(1, ii-1) = me.topLoadArray(1, ii); + end + + me.topPosArray(1, me.lengthRequired(1)) = -1.0 * polishedRodPosition / 12.0; + + if polishedRodPosition > lastPolishedRodPosition + me.topLoadArray(1, me.lengthRequired(1)) = loadMult * ... + (polishedRodLoad - me.rodWeightFluidTotal) - me.stuffingBoxFriction; + elseif polishedRodPosition < lastPolishedRodPosition + me.topLoadArray(1, me.lengthRequired(1)) = loadMult * ... + (polishedRodLoad - me.rodWeightFluidTotal) + me.stuffingBoxFriction; + else + me.topLoadArray(1, me.lengthRequired(1)) = loadMult * ... + (polishedRodLoad - me.rodWeightFluidTotal); + end + + tap=1; + while tap <= tapersAllowed + me.count(tap) = me.count(tap) + 1; + if me.count(tap) >= me.lengthRequired(tap) + if tap+1 <= me.numTapers + % working our way down to the bottom of the well + for pti = 2:me.lengthRequired(tap+1)+1 + me.topPosArray(tap+1, pti-1) = me.topPosArray(tap+1, pti); + me.topLoadArray(tap+1, pti-1) = me.topLoadArray(tap+1, pti); + end + pumpPosition = me.position(tap); + pumpLoad = me.load(tap); + status = 0; + me.topPosArray(tap+1, me.lengthRequired(tap+1)) = pumpPosition; + me.topLoadArray(tap+1, me.lengthRequired(tap+1)) = pumpLoad; + else + pumpPosition = -12.0 * me.position(tap); % + me.stretch(me.numTapers); + me.currentDownholePosition = pumpPosition; + pumpLoad = me.load(tap) + me.force(me.numTapers); + me.currentDownholeLoad = pumpLoad; + status = 1; + end + me.count(tap) = me.count(tap) - 1; + tapersAllowed = tapersAllowed + 1; + if tapersAllowed > me.numTapers + tapersAllowed = me.numTapers; + end + end + tap = tap + 1; + end + end +
Check if the stroke has ended
function directionChanged = checkEndOfStroke(me, numConsecutivePoints) + directionChanged = 0; + tempDirection = me.DIRECTION_UNKNOWN; + startDirection = me.DIRECTION_UNKNOWN; + consecutivePointsAllSameDirection = 1; + + if me.positionSensor.history(1) > me.positionSensor.history(2) + tempDirection = me.DIRECTION_UP; + startDirection = me.DIRECTION_UP; + elseif me.positionSensor.history(1) < me.positionSensor.history(2) + tempDirection = me.DIRECTION_DOWN; + startDirection = me.DIRECTION_DOWN; + end + + if startDirection == me.DIRECTION_UP + + for i = 1:numConsecutivePoints + if me.positionSensor.history(i) <= me.positionSensor.history(i+1) + consecutivePointsAllSameDirection = 0; + end + end + + if consecutivePointsAllSameDirection == 1 + tempDirection = me.DIRECTION_UP; + end + + elseif startDirection == me.DIRECTION_DOWN + + for i = 1:numConsecutivePoints + if me.positionSensor.history(i) >= me.positionSensor.history(i+1) + consecutivePointsAllSameDirection = 0; + end + end + + if consecutivePointsAllSameDirection == 1 + tempDirection = me.DIRECTION_DOWN; + end + end + + if tempDirection ~= me.lastDirection + if (tempDirection == me.DIRECTION_UP) && ... + (me.currentCard.numPointsUsed > 1) + directionChanged = 1; + end + me.lastDirection = tempDirection; + end + end +
Evaluation function for simulated position and load points
function evalSim(me, s_pos, s_load) + me.positionSensor.read(s_pos); + me.loadSensor.read(s_load); + + [d_pos, d_load, d_status] = me.calc(me.positionSensor.lastValue, ... + me.positionSensor.history(2), me.loadSensor.lastValue); + if d_status == 1 + me.currentCard.push(me.positionSensor.lastValue, ... + me.loadSensor.lastValue, d_pos, d_load); + end + + if me.checkEndOfStroke(5) == 1 + me.endOfStroke() + end + end +
End of Stroke actions
function endOfStroke(me) + me.currentCard.calcStrokeData(100, me.fluidGradient, me.rodDepthTotal, ... + me.tubingAnchorDepth, me.tubingCrossSectionalArea, ... + me.pumpArea, me.frictionEstimate, ... + me.structuralRating, me.kFactor, me.waterBblRatio, me.oilBblRatio, ... + me.gasMcfRatio); + me.plotCards() + + me.lastCard = me.currentCard(); + me.currentCard = Card(me.lastCard.strokeNumber + 1); + end +
Draws Cards
Draws current card and shows last card
function plotCards(me) + ax1 = subplot(2,2,1); + plot(me.lastCard.surfacePosition(1:me.lastCard.numPointsUsed), ... + me.lastCard.surfaceLoad(1:me.lastCard.numPointsUsed)) + title('Last Surface Card') + + ax3 = subplot(2,2,3); + plot(me.lastCard.downholePosition(1:me.lastCard.numPointsUsed), ... + me.lastCard.downholeLoad(1:me.lastCard.numPointsUsed), 'r') + title('Last Downhole Card') + + ax2 = subplot(2,2,2); + plot(me.currentCard.surfacePosition(1:me.currentCard.numPointsUsed), ... + me.currentCard.surfaceLoad(1:me.currentCard.numPointsUsed)) + title('Current Surface Card') + + ax4 = subplot(2,2,4); + plot(me.currentCard.downholePosition(1:me.currentCard.numPointsUsed), ... + me.currentCard.downholeLoad(1:me.currentCard.numPointsUsed), 'r') + title('Current Downhole Card') + + linkaxes([ax1, ax2, ax3, ax4], 'x') + + end +
end
+
+ methods(Static)
+Determines Rod Weight Per Foot given Young's Modulus and Diameter
function wtPerFt = lookupRodWeightPerFoot(i_ym, i_diam) + YM_STEEL=30.5; + YM_FIBERGLASS=7.2; + wtPerFt = -1; + if (i_ym == YM_STEEL) + if (i_diam <= 2 && i_diam > 1.75) + wtPerFt = 10.7; + elseif (i_diam <= 1.75 && i_diam > 1.65) + wtPerFt = 8.2; + elseif (i_diam <= 1.65 && i_diam > 1.5) + wtPerFt = 7; + elseif (i_diam <= 1.5 && i_diam > 1.375) + wtPerFt = 6; + elseif (i_diam <= 1.375 && i_diam > 1.125) + wtPerFt = 5; + elseif (i_diam <= 1.125 && i_diam > 1) + wtPerFt = 3.676; + elseif (i_diam <= 1 && i_diam > 0.875) + wtPerFt = 2.904; + elseif (i_diam <= 0.875 && i_diam > 0.75) + wtPerFt = 2.224; + elseif (i_diam <= 0.75 && i_diam > 0.625) + wtPerFt = 1.634; + elseif (i_diam <= 0.625 && i_diam > 0.5) + wtPerFt = 1.13; + elseif (i_diam <= 0.5) + wtPerFt = 0.72; + else + wtPerFt = 0; + end + elseif (i_ym == YM_FIBERGLASS) + if (i_diam <= 1.25 && i_diam > 1.125) + wtPerFt = 1.2879; + elseif (i_diam <= 1.125 && i_diam > 1) + wtPerFt = 1.09; + elseif (i_diam <= 1 && i_diam > 0.875) + wtPerFt = 0.8188; + elseif (i_diam <= 0.875 && i_diam > 0.75) + wtPerFt = 0.6108; + elseif (i_diam <= 0.75) + wtPerFt = 0.484; + else + wtPerFt = 0; + end + end + end +
end +end +