Files
POC-Matlab/wellSetup.m
2017-09-07 17:38:29 -05:00

208 lines
6.6 KiB
Matlab

%% Well Setup Function
function wellParams = wellSetup(...
dt, ...
rodLengths, ...
rodDiameters, ...
rodMaterials, ...
rodDampingFactors, ...
tubingHeadPressure, ...
stuffingBoxFriction, ...
fluidGradient, ...
pumpDiameter,...
tubingOuterDiameter, ...
tubingInnerDiameter, ...
tubingAnchorDepth, ...
structuralRating)
% FUNCTION PARAMETERS
% dt
% rodLengths
% rodDiameters
% rodMaterials
% tubingHeadPressure
% pumpDiameter
% tubingOuterDiameter
% tubingInnerDiameter
wellParams = struct();
wellParams.dt = dt;
wellParams.rodLengths = rodLengths;
wellParams.rodDiameters = rodDiameters;
wellParams.rodMaterials = rodMaterials;
wellParams.rodDampingFactors = rodDampingFactors;
wellParams.tubingHeadPressure = tubingHeadPressure;
wellParams.stuffingBoxFriction = stuffingBoxFriction;
wellParams.fluidGradient = fluidGradient;
wellParams.pumpDiameter = pumpDiameter;
wellParams.tubingOuterDiameter = tubingOuterDiameter;
wellParams.tubingInnerDiameter = tubingInnerDiameter;
wellParams.tubingAnchorDepth = tubingAnchorDepth;
wellParams.structuralRating = structuralRating;
numTapers = length(rodLengths);
wellParams.numTapers = numTapers;
wellParams.tubingCrossSectionalArea = (pi / 4) * ...
(power(tubingOuterDiameter, 2) - ...
power(tubingInnerDiameter,2));
wellParams.pumpArea = power(pumpDiameter, 2) * pi;
% Calculated Taper Parameters
area = zeros(1, numTapers+1);
a = zeros(1, numTapers);
pressure = zeros(1, numTapers);
rodYMs = zeros(1, numTapers);
rodWeightPerFoots = zeros(1, numTapers);
force = zeros(1, numTapers);
alpha = zeros(1, numTapers);
stretch = zeros(1, numTapers);
xOverA = zeros(1, numTapers);
lagIndex = zeros(1, numTapers, 'uint16');
factorArray = zeros(1, numTapers);
centerPoint = zeros(1, numTapers, 'uint16');
lengthRequired = zeros(1, numTapers, 'uint16');
% Calculated Rod String Totals
rodDepthTotal = 0.0;
buoyantForceTotal = 0.0;
rodWeightAirTotal = 0.0;
for i = 1:numTapers
area(i) = pi / 4 * power(rodDiameters(i), 2);
if lower(rodMaterials(i)) == "fiberglass"
rodYMs(i) = 7.2;
else
rodYMs(i) = 30.5;
end
rodWeightPerFoots(i) = lookupRodWeightPerFoot(rodYMs(i), ...
rodDiameters(i));
end
wellParams.area = area;
wellParams.rodYMs = rodYMs;
wellParams.rodWeightPerFoots = rodWeightPerFoots;
for i = 1:numTapers
a(i) = 1000.0 * sqrt(32.2 * rodYMs(i) * area(i) / rodWeightPerFoots(i));
rodDepthTotal = rodDepthTotal + rodLengths(i);
if i > 1
pressure(i) = pressure(i - 1) + fluidGradient * rodLengths(i);
else
pressure(i) = tubingHeadPressure + fluidGradient * rodLengths(i);
end
buoyantForceTotal = buoyantForceTotal + pressure(i) * (area(i+1) - area(i));
rodWeightAirTotal = rodWeightAirTotal + rodWeightPerFoots(i) * rodLengths(i);
rodWeightFluidTotal = rodWeightAirTotal + buoyantForceTotal;
end
wellParams.a = a;
wellParams.rodDepthTotal = rodDepthTotal;
wellParams.pressure = pressure;
wellParams.buoyantForceTotal = buoyantForceTotal;
wellParams.rodWeightAirTotal = rodWeightAirTotal;
wellParams.rodWeightFluidTotal = rodWeightFluidTotal;
for j = 1:numTapers
weightData = 0.0;
annularForceData = 0.0;
for i = j+1:numTapers
weightData = weightData + rodWeightPerFoots(i) * rodLengths(i);
end
for i = j:numTapers-1
annularForceData = annularForceData + pressure(i) * (area(i) - area(i+1));
end
force(j) = (-area(numTapers) * pressure(numTapers)) + weightData - annularForceData;
alpha(j) = (force(j) + rodWeightPerFoots(j) * rodLengths(j)) / (rodYMs(j) * power(10, 6) * area(j));
if j > 1
stretch(j) = stretch(j - 1) + alpha(j) * rodLengths(j) - ...
(rodWeightPerFoots(j) * power(rodLengths(j), 2.0)) / ...
(2.0 * rodYMs(j) * power(10, 6) * area(j));
else
stretch(j) = 0.0 + alpha(j) * rodLengths(j) - ...
(rodWeightPerFoots(j) * power(rodLengths(j), 2.0)) / ...
(2.0 * rodYMs(j) * power(10, 6) * area(j));
end
end
wellParams.force = force;
wellParams.alpha = alpha;
wellParams.stretch = stretch;
for i = 1:numTapers
xOverA(i) = rodLengths(i) / a(i);
lagIndex(i) = floor(rodLengths(i) / (a(i) * dt));
factorArray(i) = (xOverA(i) - double(lagIndex(i)) * dt) / dt;
centerPoint(i) = lagIndex(i) + 2;
lengthRequired(i) = 2 * (lagIndex(i) + 1) + 1;
end
wellParams.xOverA = xOverA;
wellParams.lagIndex = lagIndex;
wellParams.factorArray = factorArray;
wellParams.centerPoint = centerPoint;
wellParams.lengthRequired = lengthRequired;
frictionEstimate = rodDepthTotal * 0.10;
wellParams.frictionEstimate = frictionEstimate;
end
%% 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