Contents

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

Published with MATLAB® R2017a