From df1d2ec893100716d5a5d95a3c3df08f6d3f7186 Mon Sep 17 00:00:00 2001 From: Patrick McDonagh Date: Thu, 27 Jul 2017 20:51:22 -0500 Subject: [PATCH] Able to run simulations of cards --- AnalogInput.m | 12 ++-- AnalogInputSim.m | 50 ++++++++++++++ Card.m | 55 +++++++++------ MuxSetupSim.m | 21 ++++++ Well.m | 176 ++++++++++++++++++++++++++++++++++++++++++----- test.m | 18 +++-- 6 files changed, 284 insertions(+), 48 deletions(-) create mode 100644 AnalogInputSim.m create mode 100644 MuxSetupSim.m diff --git a/AnalogInput.m b/AnalogInput.m index f863a0d..ccf8596 100644 --- a/AnalogInput.m +++ b/AnalogInput.m @@ -13,6 +13,7 @@ classdef AnalogInput < handle properties(Access=private) m,b; + histi; end methods @@ -23,17 +24,23 @@ classdef AnalogInput < handle obj.rawMax = double(rawMax); obj.euMin = double(euMin); obj.euMax = double(euMax); + obj.histi = 0; obj.m = (obj.euMax - obj.euMin) / (obj.rawMax - obj.rawMin); obj.b = obj.euMax - obj.m * (obj.rawMax); end function value = setValue(obj, inValue) + obj.histi = obj.histi + 1; obj.rawValue = inValue; value = obj.m * inValue + obj.b; obj.lastValue = value; obj.lastStored = now; - obj.history = [value, obj.history(1:end-1)]; + + % Store value in history array + histTemp = obj.history(1:end-1); + obj.history(2:end) = histTemp; + obj.history(1) = value; end function value = read(obj) @@ -47,15 +54,12 @@ classdef AnalogInput < handle end if (obj.badReads > 10) - "Bad Reads" - obj.badReads pause(0.010); end end function value = readSim(obj, simRaw) value = obj.m * simRaw + obj.b; - obj.lastValue = pv; obj.lastValue = value; obj.lastStored = now; obj.history = [value, obj.history(1:end-1)]; diff --git a/AnalogInputSim.m b/AnalogInputSim.m new file mode 100644 index 0000000..3763f7f --- /dev/null +++ b/AnalogInputSim.m @@ -0,0 +1,50 @@ +classdef AnalogInputSim < handle + properties + mux; + channel; + rawValue; + lastValue; + lastStored=0; + rawMax,rawMin,euMax,euMin; + history=zeros(1, 100, 'double'); + badReads=0; + + end + + properties(Access=private) + m,b; + end + + methods + function obj = AnalogInputSim(mux, channel, rawMin, rawMax, euMin, euMax) + obj.mux = mux; + obj.channel = channel; + obj.rawMin = double(rawMin); + obj.rawMax = double(rawMax); + obj.euMin = double(euMin); + obj.euMax = double(euMax); + + obj.m = (obj.euMax - obj.euMin) / (obj.rawMax - obj.rawMin); + obj.b = obj.euMax - obj.m * (obj.rawMax); + end + + function value = setValue(obj, inValue) + obj.rawValue = inValue; + value = obj.m * inValue + obj.b; + obj.lastValue = value; + obj.lastStored = now; + obj.history = [value, obj.history(1:end-1)]; + end + + function value = read(obj, simRaw) + value = obj.m * simRaw + obj.b; + obj.lastValue = value; + obj.lastStored = now; + + % Store value in history array + histTemp = obj.history(1:end-1); + obj.history(2:end) = histTemp; + obj.history(1) = value; + end + end +end \ No newline at end of file diff --git a/Card.m b/Card.m index c4e60c9..4ee0736 100644 --- a/Card.m +++ b/Card.m @@ -65,7 +65,7 @@ classdef Card < handle anchorDepth, tubingCSA, pumpArea, frictionEstimate, ... structuralRating, kFactor, waterBBLRatio, oilBBLRatio, ... gasMCFRatio) - calculateSPM(); + obj.calculateSPM(); obj.surfacePositionMax = obj.positionMax(obj.surfacePosition, obj.surfaceLoad, obj.numPointsUsed); obj.surfaceLoadMax = obj.loadMax(obj.surfacePosition, obj.surfaceLoad, obj.numPointsUsed); obj.surfacePositionMin = obj.positionMin(obj.surfacePosition, obj.surfaceLoad, obj.numPointsUsed); @@ -96,21 +96,21 @@ classdef Card < handle dhLoadAtTargetTop = 0.0; dhLoadAtTargetBottom = 0.0; - for j = 1:obj.numPointsUsed + for j = 1:obj.numPointsUsed-1 if (obj.downholePosition(j) <= dhPosTarget && obj.downholePosition(j+1) > dhPosTarget) - dhLoadAtTargetTop = lineResolve(obj.downholePosition(j), obj.downholePosition(j+1), obj.downholeLoad(j), obj.downholeLoad(j+1), dhPosTarget); + dhLoadAtTargetTop = obj.lineResolve(obj.downholePosition(j), obj.downholePosition(j+1), obj.downholeLoad(j), obj.downholeLoad(j+1), dhPosTarget); end if (obj.downholePosition(j) > dhPosTarget && obj.downholePosition(j+1) >= dhPosTarget) - dhLoadAtTargetBottom = lineResolve(obj.downholePosition(j), obj.downholePosition(j+1), obj.downholeLoad(j), obj.downholeLoad(j+1), dhPosTarget); + dhLoadAtTargetBottom = obj.lineResolve(obj.downholePosition(j), obj.downholePosition(j+1), obj.downholeLoad(j), obj.downholeLoad(j+1), dhPosTarget); end if (obj.surfacePosition(j) <= suPosTarget && obj.surfacePosition(j+1) > suPosTarget) - suLoadAtTargetTop = lineResolve(obj.surfacePosition(j), obj.surfacePosition(j+1), obj.surfaceLoad(j), obj.surfaceLoad(j+1), suPosTarget); + suLoadAtTargetTop = obj.lineResolve(obj.surfacePosition(j), obj.surfacePosition(j+1), obj.surfaceLoad(j), obj.surfaceLoad(j+1), suPosTarget); end if (obj.surfacePosition(j) > suPosTarget && obj.surfacePosition(j+1) >= suPosTarget) - suLoadAtTargetBottom = lineResolve(obj.surfacePosition(j), obj.surfacePosition(j+1), obj.surfaceLoad(j), obj.surfaceLoad(j+1), suPosTarget); + suLoadAtTargetBottom = obj.lineResolve(obj.surfacePosition(j), obj.surfacePosition(j+1), obj.surfaceLoad(j), obj.surfaceLoad(j+1), suPosTarget); end end @@ -135,11 +135,11 @@ classdef Card < handle obj.downholeNetStrokeLength = obj.bottomCorner.position - obj.downholePositionMin.position; obj.fillageCalculated = (obj.downholeNetStrokeLength / obj.downholeAdjustedGrossStrokeLength) * 100.0; obj.fillageEstimated =(obj.downholeNetStrokeLength / obj.downholeGrossStrokeLength) * 100.0; - obj.fluidBBLMoved = obj.downholeNetStrokeLength * pumpArea * 0.00010307; - obj.fluidBBLMovedAdjusted = obj.fluidBBLMoved * kFactor; - obj.oilBBLMoved = obj.fluidBBLMoved * oilBBLRatio; - obj.waterBBLMoved = obj.fluidBBLMoved * waterBBLRatio; - obj.gasMCFMoved = obj.fluidBBLMoved * gasMCFRatio; + obj.fluidBblMoved = obj.downholeNetStrokeLength * pumpArea * 0.00010307; + obj.fluidBblMovedAdjusted = obj.fluidBblMoved * kFactor; + obj.oilBblMoved = obj.fluidBblMoved * oilBBLRatio; + obj.waterBblMoved = obj.fluidBblMoved * waterBBLRatio; + obj.gasMcfMoved = obj.fluidBblMoved * gasMCFRatio; if (obj.fillageEstimated > 100) @@ -169,7 +169,7 @@ classdef Card < handle obj.strokeSpeed = 60000.0 / (now - obj.strokeStartTime); end - function dist = distanceToLine(pos, load) + function dist = distanceToLine(obj, pos, load) % Finds the distance between a point and a line between the % Max Position at Max Load and Min Position at Min Load x1 = obj.downholePositionMin.position; @@ -180,6 +180,27 @@ classdef Card < handle dist = abs((y2-y1)* pos - (x2-x1)* load + x2*y1 - y2*x1) / sqrt(power(y2-y1, 2) + power(x2-x1,2)); end + function plot(obj) + ax1 = subplot(2,1,1); + plot(obj.surfacePosition(1:obj.numPointsUsed), obj.surfaceLoad(1:obj.numPointsUsed)) + title('Surface Card') + + ax2 = subplot(2,1,2); + plot(obj.downholePosition(1:obj.numPointsUsed),obj.downholeLoad(1:obj.numPointsUsed), 'r') + title('Downhole Card') + + linkaxes([ax1, ax2], 'x') + end + end + + methods(Static) + function yTest = lineResolve(x1, x2, y1, y2, xTest) + % Uses the standard line equation to find the y value given x + line_m = (y2 - y1) / (x2 - x1); + line_b = y1 - line_m * x1; + yTest = line_m * xTest + line_b; + end + function foundPoint = positionMax(positionArr, loadArr, arrSize) maxPos = positionArr(1); loadAtMaxP = -inf; @@ -228,16 +249,6 @@ classdef Card < handle end foundPoint = LPPair(posAtMinL, minLoad); end - - end - - methods(Static) - function yTest = lineResolve(x1, x2, y1, y2, xTest) - % Uses the standard line equation to find the y value given x - line_m = (y2 - y1) / (x2 - x1); - line_b = y1 - line_m * x1; - yTest = line_m * xTest + line_b; - end end diff --git a/MuxSetupSim.m b/MuxSetupSim.m new file mode 100644 index 0000000..7987a94 --- /dev/null +++ b/MuxSetupSim.m @@ -0,0 +1,21 @@ +classdef MuxSetupSim < handle + properties(Access=private) + mux1Pin=5; + mux2Pin=6; + mux3Pin=13; + digInPin=19; + anOutTriggerPin=23; + dev; + end + + properties + setup=[0 0 0; 1 0 0; 0 1 0; 1 1 0; 0 0 1; 1 0 1; 0 1 1; 1 1 1]; + spiDevice; + end + + methods + function obj = MuxSetupSim(dev) + obj.dev = dev; + end + end +end \ No newline at end of file diff --git a/Well.m b/Well.m index 4f8cae6..48c4b3a 100644 --- a/Well.m +++ b/Well.m @@ -3,27 +3,49 @@ classdef Well < handle % 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); + currentCard=Card(0); + lastCard=Card(0); + direction=0; + lastDirection=0; % User Inputs - dt = 0.03333333; + 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; @@ -31,7 +53,6 @@ classdef Well < handle rodDiameter; rodYM; rodWeightPerFoot; - % Calculated Taper Parameters a; @@ -47,16 +68,17 @@ classdef Well < handle centerPoint; count; - topPosArray; - topLoadArray; - 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; @@ -64,13 +86,10 @@ classdef Well < handle loadBefore3 = 0.0; loadAfter3 = 0.0; tempLoadValue = 0.0; - - - end methods - function me = Well(name) + function me = Well(name, sim) me.wellName = name; me.c = [0.8 0.8 0.8]; @@ -97,10 +116,27 @@ classdef Well < handle me.topPosArray = zeros(10, 250,'double'); me.topLoadArray = zeros(10, 250, 'double'); - me.taperSetup(); + 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 taperSetup(me) + 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 @@ -153,6 +189,8 @@ classdef Well < handle 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) @@ -174,7 +212,7 @@ classdef Well < handle me.loadAfter3 = 0.0; pumpPosition = exp(ci * lengthi / (2.0 * ai)) * (me.topPosArray(i, iAfter) + ... - factori * me.topPosArray(i, iAfter+1)); + 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) - ... @@ -271,6 +309,8 @@ classdef Well < handle function [pumpPosition, pumpLoad, status] = calc(me, polishedRodPosition, lastPolishedRodPosition,... polishedRodLoad) + me.currentSurfacePosition = polishedRodPosition; + me.currentSurfaceLoad = polishedRodLoad; pumpPosition = -1.0; pumpLoad = -1.0; status = -1; @@ -312,7 +352,9 @@ classdef Well < handle 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; @@ -324,6 +366,108 @@ classdef Well < handle 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) diff --git a/test.m b/test.m index 4e9583b..116740a 100644 --- a/test.m +++ b/test.m @@ -1,10 +1,16 @@ malletData = csvread('mallet.csv'); -mallet = Well("Mallet"); -for i = 2:size(malletData,1) - [dhPos, dhLoad, status] = mallet.calc(malletData(i,1), malletData(i-1,1), malletData(i,2)); - if status == 1 - mallet.currentCard.push(malletData(i,1), malletData(i,2), dhPos, dhLoad); +mallet = Well("Mallet", true(1)); +figure; +for loops = 1:6 + for i = 2:size(malletData,1) + mallet.evalSim(malletData(i,1), malletData(i,2)); end end -figure;plot(mallet.currentCard.downholePosition,mallet.currentCard.downholeLoad) + +% mallet.currentCard.calcStrokeData(100, mallet.fluidGradient, mallet.rodDepthTotal, ... +% mallet.tubingAnchorDepth, mallet.tubingCrossSectionalArea, ... +% mallet.pumpArea, mallet.frictionEstimate, ... +% mallet.structuralRating, mallet.kFactor, mallet.waterBblRatio, mallet.oilBblRatio, ... +% mallet.gasMcfRatio) +