clear
clc

%**************************************************************************
% MANUAL INPUT
Pr = 250        % Reservoir pressure

Pb = 240        % Bubblepoint pressure
Pwh = 30        % Wellhead pressure (Either sep or seafloor depending on "PandR")
pump = 0        % 1 = pump on 0 = pump off
Pp = 10         % Pump output pressure
PandR = 0       % 1 = add pipiline and riser(1/0)
Pplace = 300    % Pump placement =  number * DL m from foot
GORS = 100;     % Seperator GOR [Sm3/Sm3]
%**************************************************************************

% Input data list

g = 9.81;               % Gravity
Lh = 1000;              % length of horizontal well section
Lv = 2000;              % Length of vertical section
Lp = 2000;              % Length pipeleine
Lr = 300;               % Length riser
Dl = 10;                % Segment length
K = 1.3e-3;             % Perforation inflow constant
Rol = 800;              % Oil density
Rog0 = 1.5;             % Gas density at reference conditions
Myl = 2e-3;             % liquid viscosity
Nyg = 2e-5;             % Gas viscocity
Dhor = 0.15;            % Diameter horizontal section
Dver = 0.20;            % Diameter vertical section
Dpip = Dver;            % Diameter vertical section
Dris = Dver;            % Diameter vertical section
angleh = 0;
anglev = pi/2;
anglep = 0;
angler = pi/2;
sig = 30;
Pwiguess = Pr-15;
Pwi(1) = Pwiguess;
DT = 30;                   % Temperature increase /km
if PandR ==1
    Ltot = Lh+Lv+Lp+Lr;
else
    Ltot = Lh+Lv;
end




%**************************************************************************
% Pre analysis.
% Determine wich cells in the horizontal section wich contains perforations

Perfor = linspace(0,800,5);
Nper = length(Perfor);
Nhor = 1+Lh/Dl;                % Number of cells in horizontal
Nver = Lv/Dl;                  % Number of vertical cells
Npip = Lp/Dl;                  % Number of pipeline cells
Nris = Lr/Dl;                  % Number of riser cells
Ntot = Nhor+Nver;
if PandR == 1
    Ntot = Ntot+Npip+Nris;
end
Kxi = zeros(1,Ntot);

%**************************************************************************
%ITERATION
noit = 1;       %Iterations
ac = 0.0001;   % Accuracy
error(1) = 1;

while abs(error(noit)) > ac

        Pwi(1) = Pwi(1) - 0.1*error(noit);

        noit = noit + 1;

%**************************************************************************

for j = 1:Nhor
       Lpos(j) = (j-1)*Dl;
    for iper = 1:Nper
        if Lpos(j) == Perfor(iper)
            Kxi(j) = K;         % Set performation "flag": Inflow constant
        end
    end
end

for k = 1:Ntot                  % Make angle and diameter vectors
    if k <= Nhor
        angle(k) = angleh;
        D(k) = Dhor;
    elseif k > Nhor && k <= Nver
        angle(k) = anglev;
        D(k) =Dver;
    elseif k > Nver && k <= Npip
        angle(k) = anglep;
        D(k) = Dpip;
    else
        angle(k) = angler;
        D(k) = Dris;
    end
end

%**************************************************************************

Mtot = 0;                           % Mass flow rate

for i = 1:Ntot
    if Kxi(i) > 0
        if Pwi(i) > Pb
        Minn = Rol*K*(Pr-Pwi(i));   % Mass inflow through perf.
        else
        Pvogel = Pwi(i);
        Minn = vogel(Pvogel,K,Pr,Pb,Rol); % Vogel inflow
        end
        Mtot = Mtot+Minn;
    end

%*************************************************************************
% GAS FRACTION AND TEMPERATURE CALCULATION, Hasan-Kabir Correlation

%%%%%TEMPERATURE HASSAN-KABIR
% Note that values are converted to match input parameters, then converted
% back

%Input
kcem = 0.42;
ks = 26;
kesa = 1.06;
kesh = 0.7;
ko = 0.08;                                  % Btu/hr-ft-F
Cpo = 0.65;                                 % Btu/lbm-F
alphasa = 0.04;
alphash = 0.02;
Tinitial = 277.15+DT*Lv*10^(-3);            % Reservoir temperature
rto =	2.75;                               % In
rti= 2.446;                                 % In
rw = 4;                                     % In
gG1 = 0.0083;
gG2	= 0.015;
hdi = 2.45;
TVD = 2000;
depthi = 1800;
depthsa = 200;
t = 24*365*12;                              % Time [hours]
rohw = 62.4;                                % lbm/ft^3
So = Rol/1000;                              % Frac

MinnFU = Minn*2.2046;
myoFU = exp(7.296*(((So)^3)-3.095));        % cP
myo = myoFU*(((10^(-3))*2.2046)/(0.3048^(-1))); %lbm/ft-sec

tdwsa = (alphasa*t)/(rw^2);
tdwsh = (alphash*t)/(rw^2);

if tdwsa <= 1.5
    ftsa = 1.1281*sqrt(tdwsa)*(1-0.3*sqrt(tdwsa));
else
    ftsa = 0.4063+0.5*log(tdwsa)*(1+(0.6/tdwsa));
end
if tdwsh <= 1.5
    ftsh = 1.1281*sqrt(tdwsh)*(1-0.3*sqrt(tdwsh));
else
    ftsh = 0.4063+0.5*log(tdwsh)*(1+(0.6/tdwsh));
end

Nnu = 0.023*(((rohw*So*2*rti)/myo)^0.8)*((3600*Cpo*myo)/ko)^(1/3);
hf = (Nnu*2*rti)/ko;

if i < 130
Uo = ((rto/(rti*hf))+((rto*log(rto/rti))/ks)+((rto*log(rw/rto))/kcem)+(rto*ftsa)/kesa)^(-1);
else
Uo = ((rto/(rti*hf))+((rto*log(rto/rti))/ks)+((rto*log(rw/rto))/kcem)+(rto*ftsh)/kesh)^(-1);
end

A = (MinnFU*Cpo)/(pi*2*rto*Uo);

Tei(1) = (Tinitial-273.15)*9/5 + 32;
Tf(1) = (Tinitial-273.15)*9/5 + 32;
Tei(i+1) = Tei(i)-gG1*Dl*sin(angle(i));

if i < 131 %Sandstone
    Tf(i+1) = Tei(i)-gG1*Dl*sin(angle(i))+(Tf(i)-Tei(i))*exp(-Dl/A)+ gG1*sin(angle(i))*A*(1-(exp(-Dl/A)));
else % Shale
    Tf(i+1) = Tei(i)-gG2*Dl*sin(angle(i))+(Tf(i)-Tei(i))*exp(-Dl/A)+ gG2*sin(angle(i))*A*(1-(exp(-Dl/A)));
end

T(i) = (Tf(i)-32)*5/9+273.15;



%%%%%%%GAS FRACTION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%z = nonlinear(Pwi(i),T(i));
z = 1;
C = (Rog0*298.15);
Rog(i) = (Pwi(i)*C)/(z*T(i));
GOR(i) = GORS+10^-7+((GORS*(Pwh-Pwi(i))/(Pwi(1)-Pwh)));
Gfrak(i) = (1+Rol/(Rog(i)*GOR(i)))^(-1);

%*************************************************************************
    if Pwi(i) > Pb
        Rom = Rol;
        Mym = Myl;
        Ulsi(i) = Mtot/(Rol*(pi*D(i)^2/4));
        Umix = Ulsi(i);
        Epsgi(i) = 0;
        Gfrak(i) = 0;


        Reyn = Rom*Umix*Dhor/Mym;
        Frik = 0.046*Reyn^(-0.2);
        Dpdxf = (4/Dhor)*Frik*.5*Rom*Umix^2;
        Pwi(i+1) = Pwi(i) - Dpdxf*Dl*1e-5;

    else



        Mgi(i) = Gfrak(i)*Mtot;
        Mli(i) = Mtot-Mgi(i);
        Ugsi(i) = Mgi(i)/(Rog(i)*(pi*D(i)^2/4));
        Ulsi(i) = Mli(i)/(Rol*(pi*D(i)^2/4));
        Myg = Nyg*Rog(i);

        %For plot purpose (Calculatedin bb)
        Umix = Ulsi(i)+Ugsi(i);
        Epsgi(i) = Ugsi(i)/Umix;
       %************************

        Pwi(i+1) = Pwi(i) - bb(Ulsi(i),Ugsi(i),Rol,Rog,Myl,Myg,sig,D(i),angle(i),g,Pwi(i),Dl)*1e-5;

        if pump == 1 && i == Pplace
        Pwi(i+1) = Pwi(i+1) + Pp;
        end
    end

end
Pwhout = Pwi(Ntot)
error(noit) = (Pwhout-Pwh)/(Pwh);

Mtot


figure(1)
subplot(3,1,1); plot(linspace(0,(Ltot),length(Pwi)),Pwi,'k-')
title('Pressure vs. Position')
ylabel('Pressure [bar]')
xlabel('Position from toe [m]')
subplot(3,1,2); plot(linspace(0,(Ltot),length(Epsgi)),Epsgi,'k-')
title('Epsgi')
ylabel('Fraction')
xlabel('Position from toe [m]')
subplot(3,1,3); plot(linspace(0,noit,noit),error,'k-')
ylabel('Error')
xlabel('Iterations')
drawnow

end


Published with MATLAB® R2014b

Beggs and Brill Pressure drop correlation

Beggs and Brill Pressure drop correlation

function[dp] = bb(Ulsi, Ugsi, Rol, Rog, Myl, Myg, sig, D, angle, g, Pwi,Dl)

Nlv = Ulsi*(Rol/(g*sig))^0.25;
Nfr = ((Ulsi+Ugsi)^2)/(g*D);
Lambdal = Ulsi/(Ulsi+Ugsi);
Lambdag = Ugsi/(Ulsi+Ugsi);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Determination of flow regime

L1 = 316*Lambdal^0.302;
L2 = 9.52*(10^(-4))*Lambdal^(-2.4684);
L3 = 0.1*Lambdal^(-1.4516);
L4 = 0.5*Lambdal^(-6.738);

FC = 0;                 % Flow regime controller
% FC: 1 = Segregated, 2 = Transition, 3 = Intermittent, 4 = Distributed

if Lambdal < 0.01 && Nfr < L1
    FC = 1;
end
if Lambdal >= 0.01 && Nfr < L2
    FC = 1;
end
if  Lambdal >= 0.01 && L2 <= Nfr && Nfr <= L3
    FC = 2;
end
if Lambdal >= 0.01 && Lambdal < 0.4 && Nfr > L3 && Nfr <= L1
    FC = 3;
end
if Lambdal >= 0.4 && Nfr > L3 && Nfr <= L4
    FC = 3;
end
if Lambdal < 0.4 && Nfr >= L1
    FC = 4;
end
if Lambdal >= 0.4 && Nfr >= L4
    FC = 4;
end

% In case of transitional flow the following parameters are needed
A  =(L3-Nfr)/(L3-L2);
B = 1-A;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Table of constants
% Vertical
% Segregated flow
as = 0.98;
bs = 0.4846;
cs = 0.0868;
% Intermittent
ai = 0.845;
bi = 0.5351;
ci = 0.0173;
% Distributed
ad = 1.065;
bd = 0.5824;
cd = 0.0609;

% Horizontal
%Segregated upwards
ddashs = 0.011;
es = -3.768;
fs = 3.539;
gs = -1.624;
%Intermittent upwards
ddashi = 2.96;
ei = 0.305;
fi = -0.4473;
gi = 0.0978;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Liquid fraction horizontal

Hl0s = (as*Lambdal^bs)/(Nfr^cs);
Hl0i = (ai*Lambdal^bi)/(Nfr^ci);
Hl0d = (ad*Lambdal^bd)/(Nfr^cd);
Hl0t = A*Hl0s + B*Hl0i;

if FC == 1
    ddash = ddashs;
    e = es;
    f = fs;
    g = gs;
    Hl = Hl0s;
elseif FC == 2
    Hl = Hl0t;
elseif FC == 3
    ddash = ddashi;
    e = ei;
    f = fi;
    g = gi;
    Hl = Hl0i;
elseif FC == 4
    Hl = Hl0d;
else
    disp('Error in Beggs and Brill - no flow regime calculation')
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Liquid fraction inclined
if angle ~= 0 && FC ~= 4 && FC ~= 2;
    C = (1-Lambdal)*log(ddash*(Lambdal^e)*(Nlv^f)*(Nfr^g));
    Fangle = sin(1.8*angle)-(1/3)*(sin(1.8*angle))^3;
    yps = 1+ C*Fangle;
    Hl = Hl*yps;
elseif angle ~= 0 && FC == 2 % Transition inclined
    C1 = (1-Lambdal)*log(ddashs*(Lambdal^es)*(Nlv^fs)*(Nfr^gs));
    C2 = (1-Lambdal)*log(ddashi*(Lambdal^ei)*(Nlv^fi)*(Nfr^gi));
    Fangle = sin(1.8*angle)-(1/3)*(sin(1.8*angle))^3;
    yps1 = 1+ C1*Fangle;
    yps2 = 1+ C2*Fangle;
    Hl = A*Hl0s*yps1+B*Hl0i*yps2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Frictional pressure drop
y = Lambdal/Hl^2;
if y > 1 && y <1.2
    s = log(2.2*y-1.2);
else
    s = log(y)/(-0.052+3.182*log(y)-0.8725*(log(y)^2)+0.01853*(log(y))^4);
end
Ron = Rol*Lambdal+Rog*Lambdag;
myn = Myl*Lambdal+Myg*Lambdag;
Nre = (Ron*(Ulsi+Ugsi)*D)/myn;
fn = (2*log10(Nre/(4.5223*log10(Nre)-3.8215)))^(-2);
ftp = fn*exp(s);
dpf = (ftp*Ron*(Ugsi+Ulsi)^2)/(2*D);

%Hydrostatic Pressure Drop
Ros = Rol*Hl+Rog*(1-Hl);
if angle == 0
    dpx = 0;
else
    dpx = Ros*g*sin(angle);
end
%Kinetic Pressure Drop
Ek = (Ros*(Ulsi+Ugsi)*Ugsi)/(Pwi*1e5);

%Total pressure drop
dp = ((dpx+dpf)*Dl)/(1-Ek);

end


Published with MATLAB® R2014b




      
   Vogel inflow
Vogel inflow
function Minn = vogel(Pvogel,K,Pr,Pb,Rol)

%Minn at Pb for 1. perf = 10.3996 kg/s from Darcy

qmax = K*(Pr-Pb+(Pb/1.8));
qb = K*(Pr-Pb);
qo = (qmax-qb)*(1-0.2*(Pvogel/Pb)-0.8*(Pvogel/Pb)^2)+qb;
Minn = qo*Rol;

end


Published with MATLAB® R2014b




      
   nonlinear
Contents
COMPRESSIBILITY FACTOR
function[z] = nonlinear(Pressure,Temperature)

P = Pressure*14.5;
T = Temperature*(9/5);
Sg = 1.22;
%*****
J = 0.070729*Sg -3.2191*Sg^2;
K = 17.438*Sg-3.2191*Sg^2;
Tpc = (K^2)/J;
Ppc = Tpc/J;
Tred = T/Tpc; % Make sure after T calc
Pred = P/Ppc;

% Constants
A1 = 0.3265;
A2 = -1.07;
A3 = -0.5339;
A4 = 0.01569;
A5 = -0.05165;
A6 = 0.5475;
A7 = -0.7361;
A8 = 0.1844;
A9 = 0.1056;
A10 = 0.6134;
A11 = 0.7210;
%*************::

ac = 10^-4;
a = 100;
M = 1;
noit = 0;


while M > ac

    noit = noit + 1;
    z = linspace(0.1,2,a);

        for i = 1:length(z)

            rohred(i) = (0.27*Pred)/(z(i)*Tred);

            f1(i) = (A1+(A2/Tred)+(A3/(Tred^3))+(A4/(Tred^4))+(A5/(Tred^5)))*rohred(i);
            f2(i) = (A6+(A7/Tred)+(A8/(Tred^2)))*rohred(i)^2;
            f3(i) = A9*((A7/Tred)+(A8/(Tred^2)))*rohred(i)^5;
            f4(i) = A10*(1+A11*rohred(i)^2)*((rohred(i)^2)/(Tred^3))*exp(-A11*rohred(i)^2);

            f(i) = 1+f1(i)+f2(i)-f3(i)+f4(i)-z(i);
            g(i)= abs(f(i));

        end

        if min(f) > 0 % Error in calculation, set comp = 1 and break loop
                    z = 1;
                    break
        end

    [M,I] = min(g);
    z = z(I);

    a = a*10;

end
end


Published with MATLAB® R2014b