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 wellName; % Instrumentation device; mux; positionSensor; loadSensor; % Current Values currentSurfacePosition,currentSurfaceLoad; currentDownholePosition,currentDownholeLoad; currentCard=Card(0); lastCard=Card(0); direction=0; lastDirection=0; % User Inputs 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 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 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 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 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 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 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 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 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 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) 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