this code was sabotaged please correct it

2 views (last 30 days)
noobie
noobie on 27 May 2023
Commented: Walter Roberson on 5 Jun 2023
Ran in:
%%%%% Spacer 2 changed by tail cement (vsp2 Vsa2 Lsp2 Lspf2 ...) and cement (Lead) (Lcp Vcp Lcpf...)
% Conventional cemeting
% clear all
close all
clc
%Tuning
Lfactor = 1;
Ffactor =1.2;
%%%%%%%
in2m = 0.0254; % multiplier from in to m
%%%%%%%%%%%%% surface line
L_surf = 30; % length of surface line [m]
L_surf = L_surf*Lfactor;
Surfd_i = 1.69; % Surface pipe diameter [in]
Asurf = pi*(Surfd_i*in2m/2)^2; % inner area of surface pipe [m^2]
V_tsurf = Asurf*L_surf; % total volume surface pipe [m^3]
%%%%%%%%%%%%%
L_SI = 2311; % length of well, [m]
rho_mud = 1340; % mud density bentonite slurry, [kg/m^3]
rho_cem = 1500; % Portland cement, [kg/m^3]
rho_spacer = 1400; % spacer density, [kg/m^3]
rho_spacer2 = 1900;
Cd_o = 13.375; % casing outer diameter, [in]
Sd_i = 12.415; % stinger inner diameter, [in]
Od = (16*(2311-398)+(17.755*398))/2311; % Openhole diameter, [in]
PV_mud = 20; % plastic viscosity mud, [cP]
PV_spacer = 26; % plastic viscosity spacer, [cP]
PV_cem = 24; % plastic viscosity cement, [cP]
PV_spacer2 =84;
YP_SI_mud = 15; % yield point mud, [Pa]
YP_SI_spacer= 14; % yield point spacer, [Pa]
YP_SI_cem = 15; % yield point cement, [Pa]
YP_SI_spacer2 =17;
tau_y_mud = 0; % yield stress mud, [lb/100ft^2]
tau_y_spacer= 0; % yield stress spacer, [lb/100ft^2]
tau_y_cem = 0; % yield stress cement, [lb/100ft^2]
tau_y_spacer2 = 0;
g = 9.81; % gravity constant, [m/s^2]
in2m = 0.0254; % multiplier from in to m
lpm2m3s = 1/60000; % multiplier from lpm to m3/s
bar2Pa = 1e5; % multiplier from bar to Pa
B_mud = 1.5*10^4 * bar2Pa; % isothermal bulk modulus mud [Pa]
B_spacer = 1.5*10^4*bar2Pa; % isothermal bulk modulus spacer [Pa]
B_cem = 40*10^4 * bar2Pa; % isothermal bulk modulus cement [Pa]
B_spacer2 = 20*10^4 * bar2Pa;
Ap = pi*(Sd_i*in2m/2)^2; % inner area of stinger [m^2]
Aa = pi*((Od*in2m/2)^2-(Cd_o*in2m/2)^2); % inner area of annulus [m^2]
V_tp = Ap*L_SI; % total volume pipe [m^3]
V_ta = Aa*L_SI; % total volume annulus [m^3]
Vst1 = 14; % pre-flush-spacer [m^3]
Vst2 = 75.56; % total cement volume [m^3]
%
t_m = 0; % time pumping mud before cement
dt = 0.05; % time step [s]
T = 370*60; % total duration of operation [s]
t = dt:dt:T;
N = length(t);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
qp_mud = 1000*lpm2m3s; % displacing mud rate [m^3/s]
qp_spacer1 = 1000*lpm2m3s;
qp_cem = 900*lpm2m3s; % pumping cement rate [m^3/s]
qp_tail = 900*lpm2m3s;
qp_mud_dis1 = 1200*lpm2m3s; % displacement mud rate [m^3/s]
qp_mud_dis2 = 300*lpm2m3s; % displacement mud rate [m^3/s]
qp_mud_dis3 = 1000*lpm2m3s; % displacement mud rate [m^3/s]
qp_mud_dis4 = 800*lpm2m3s; % displacement mud rate [m^3/s]
qp_mud_dis5 = 300*lpm2m3s;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
P0 = 0; % gauge Pressure pressure = 0 [bar]
Lmp0 = L_SI; % initial length of mud in stinger [m]
Lma0 = Lmp0; % initial length of mud in annulus [m]
Lcem_an = 1018; % length of cement in annulus [m]
Vct = 56; % Total volume of cement to be pumped [m3]
Vanneau = 1.94 ; %Volome anneau [m3]
%% Pre-allocation of times series: a=annulus, p=pipe, c=cement, s=spacer, m=mud, f=front
zz= zeros(1,N);
Fsurf = zeros(1,N);
qp = zeros(1,N);
Ga = zeros(1,N);
Gp = zeros(1,N);
Mp = zeros(1,N);
Ma = zeros(1,N);
qb = zeros(1,N);
qout = zeros(1,N);
qavg_pipe = zeros(1,N);
Fp = zeros(1,N);
Fa = zeros(1,N);
Pp = zeros(1,N);
Ppump = zeros(1,N);
Pdha = zeros(1,N);
deltaP = zeros(1,N);
Lcp = zeros(1,N);
Lcpf = zeros(1,N);
Lca = zeros(1,N);
Lma = zeros(1,N);
Lmp = zeros(1,N);
Lsp = zeros(1,N);
Lspf = zeros(1,N);
Lsp1 = zeros(1,N);
Lsp2 = zeros(1,N);
Lsa = zeros(1,N);
Lspf2 = zeros(1,N);
Lsa2 = zeros(1,N);
hd = zeros(1,N);
Vcp = zeros(1,N);
Vcp_c = zeros(1,N);
Vmp1 = zeros(1,N);
Vmp2 = zeros(1,N);
Vmp_c1 = zeros(1,N);
Vmp_c2 = zeros(1,N);
Vma_c = zeros(1,N);
Vca_c = zeros(1,N);
Vca = zeros(1,N);
Vma = zeros(1,N);
Vsp1 = zeros(1,N);
Vsp_c1 = zeros(1,N);
Vsp2 = zeros(1,N);
Vsp_c2 = zeros(1,N);
Vsa2 = zeros(1,N);
Vsa_c2 = zeros(1,N);
Vsa = zeros(1,N);
Vsa_c = zeros(1,N);
statec = zeros(1,N);
%% Initial conditions
qp(1) = qp_mud;
qb(1) = qp_mud;
qout(1) = qp_mud;
qavg_pipe(1) = qp_mud;
Ga(1) = rho_mud*g*L_SI;
Gp(1) = rho_mud*g*L_SI;
Mp(1) = rho_mud*L_SI/((pi*(Sd_i^2)/4)*in2m^2);
Ma(1) = rho_mud*L_SI/((pi*((Od^2)-(Cd_o^2))/4)*in2m^2);
Fp(1) = Ffactor*unified_friction_pipe(Sd_i,qp_mud,L_SI,PV_mud,YP_SI_mud,tau_y_mud,rho_mud);
Unrecognized function or variable 'unified_friction_pipe'.
Fa(1) = Ffactor*unified_friction_ann(Cd_o,Od,qp_mud,L_SI,PV_mud,YP_SI_mud,tau_y_mud,rho_mud);
Fsurf(1) = Ffactor*unified_friction_pipe(Surfd_i,qp_mud,L_surf,PV_mud,YP_SI_mud,tau_y_mud,rho_mud);
Pdha(1) = P0 + Fa(1) + Ga(1);
Pp(1) = Pdha(1) + Fp(1) - Gp(1) ;
Ppump(1) = Pp(1) + Fsurf(1);
r0 = Pdha(1);
Lma(1) = Lma0;
Lmp(1) = Lmp0;
Lmp1(1) = Lmp0;
Vmp1(1) = V_tp;
Vma(1) = V_ta;
%% Flow
tramp = 20; % ramp up time rig pump
headindp = true; % true/false cement head in pipe
Tmin(1)=0;
for k = 2:N
k;
Tmin(k)=t(k)/60;
if t(k) < 14*60 % ramping up pump pressure 10min
qp(k) = qp_mud; %spacer1 1000 lpm
elseif t(k) < 24*60
qp(k) = max(0, qp(k-1)-dt*qp_mud/tramp); %shutdown 10min
elseif t(k) < 170*60
qp(k) = max(qp_tail, qp(k-1)-dt*qp_mud/tramp); %tail 900 lpm 84.5 min
elseif t(k) <182*60
qp(k) = max(0, qp(k-1)-dt*qp_mud/tramp); %shutdown 10min
%Ffactor =1.85;
elseif t(k) < 219*60
qp(k) = min(qp_mud_dis1, qp(k-1)+dt*qp_mud/tramp); % disp1
elseif t(k) < 229*60
qp(k) = max(qp_mud_dis2, qp(k-1)-dt*(qp_mud/tramp)); % disp 2
elseif t(k) < 347*60
qp(k) = min(qp_mud_dis3, qp(k-1)+dt*(qp_mud/tramp)); % disp 3
elseif t(k) < 360*60
qp(k) = min(qp_mud_dis4, qp(k-1)+dt*(qp_mud/tramp)); % disp 4
elseif t(k) < 370*60
qp(k) = max(qp_mud_dis5, qp(k-1)-dt*(qp_mud/tramp)); % disp 5
else
qp(k) = max(0, qp(k-1)-dt*qp_mud/tramp);
end
%State machine
if statec(k-1) == 0
if t(k) <= 0.00 % pumping mud
statec(k) = 0;
else
statec(k) = 1;
end
elseif statec(k-1) == 1
if Vsp1(k-1) < Vst1 %pumping spacer
statec(k) = 1;
else
statec(k) = 2;
end
elseif statec(k-1) == 2
if Vcp(k-1) < Vct %pumping cement (Lead)
statec(k) = 2;
else
statec(k) = 3;
end
elseif statec(k-1) == 3
if Vsp2(k-1) < Vst2 %pumping spacer2 (tail)
statec(k) = 3;
else
statec(k) = 4;
end
elseif statec(k-1) == 4 %pumping mud
if Lspf(k-1) < L_SI % spacer1
statec(k) = 4;
else
statec(k) = 5;
end
elseif statec(k-1) == 5
if Lcpf(k-1) < L_SI % Cement(lead) front reach bottom of well
statec(k) = 5;
else
statec(k) = 6;
end
elseif statec(k-1) == 6
if Lspf2(k-1) < L_SI % Cement(tail) front reach bottom of well
statec(k) = 6;
else
statec(k) = 7;
end
elseif statec(k-1) == 7
if Vsp2(k-1) < Vanneau % Displace with mud
statec(k) = 7;
else
statec(k) = 8;
kend = N;
end
elseif statec(k-1) == 8
statec(k) = 8; % Ramp down rig pump
end
if statec(k) >= 4 % Finished pumping cement and spacer
qp(k) = min(qp_mud, qp(k-1)+dt*qp_mud/tramp);
%
end
if statec(k) == 0
Vmp1(k) = V_tp;
Vma(k) = V_ta;
Vma_c(k) = V_ta;
Vmp_c1(k) = V_tp;
Vcp(k) = 0;
Lcp(k) = Vcp(k)/Ap;
Lcpf(k) = Lcp(k);
Pp(k) = max(0,Pp(k-1) + dt*(B_mud/V_tp)*(qp(k-1) - qb(k-1)));
Pdha(k) = max(0,Pdha(k-1) + dt*(B_mud/V_ta)*(qb(k-1) - qout(k-1)));
elseif statec(k) == 1
Vmp_c1(k) = Vmp1(k-1) - dt*(qb(k-1));
Vsp_c1(k) = Vsp1(k-1) + dt*(qp(k-1));
const_mud_pipe1 = Vmp_c1(k)/B_mud;
const_spacer_pipe1 = Vsp_c1(k)/B_spacer;
deltaP(k) = -(V_tp-hd(k-1)*Ap-Vmp_c1(k)-...
Vsp_c1(k))/(const_mud_pipe1+const_spacer_pipe1);
Vmp1(k) = Vmp_c1(k)*(1-deltaP(k)/B_mud);
Vsp1(k) = Vsp_c1(k)*(1-deltaP(k)/B_spacer);
Vma(k) = V_ta;
Vma_c(k) = V_ta;
Lsp1(k) = Vsp1(k)/Ap;
Lspf(k) = Lsp1(k) +hd(k-1);
const_mud_pipe2 = Vmp1(k)/B_mud;
const_spacer_pipe2 = Vsp1(k)/B_spacer;
Fsurf(k) = Ffactor*unified_friction_pipe(Surfd_i,qp(k-1),L_surf,PV_spacer,YP_SI_spacer,tau_y_spacer,rho_spacer);
Pp(k) = max(0,Pp(k-1) + dt*(qp(k-1)-qb(k-...
1))/(const_mud_pipe2+const_spacer_pipe2));
Pdha(k) = max(0,Pdha(k-1) + dt*(B_mud/V_ta)*(qb(k-1) - qout(k-1)));
Ppump(k) = max(0,Pp(k)+Fsurf(k));
elseif statec(k) == 2
Vmp_c1(k) = Vmp1(k-1) - dt*(qb(k-1));
Vcp_c(k) = Vcp(k-1) + dt*(qp(k-1));
Vsp_c1(k) = Vsp_c1(k-1);
const_mud_pipe1 = Vmp_c1(k)/B_mud;
const_cem_pipe1 = Vcp_c(k)/B_cem;
const_spacer_pipe1 = Vsp_c1(k)/B_spacer;
deltaP(k) = -(V_tp-hd(k-1)*Ap-Vmp_c1(k)-Vcp_c(k)-Vsp_c1(k))/...
(const_mud_pipe1+const_cem_pipe1+const_spacer_pipe1);
Vmp1(k) = Vmp_c1(k)*(1-deltaP(k)/B_mud);
Vcp(k) = Vcp_c(k)*(1-deltaP(k)/B_cem);
Vsp1(k) = Vsp_c1(k)*(1-deltaP(k)/B_spacer);
Vma(k) = V_ta;
Vma_c(k) = V_ta;
Lcp(k) = Vcp(k)/Ap;
Lcpf(k) = Lcp(k)+hd(k-1);
Lsp1(k) = Vsp1(k)/Ap;
Lspf(k) = Lcpf(k)+Lsp1(k);
const_mud_pipe2 = Vmp1(k)/B_mud;
const_cem_pipe2 = Vcp(k)/B_cem;
const_spacer_pipe2 = Vsp1(k)/B_spacer;
Fsurf(k) = Ffactor*unified_friction_pipe(Surfd_i,qp(k-1),L_surf,PV_cem,YP_SI_cem,tau_y_cem,rho_cem);
Pp(k) = max(0,Pp(k-1) + dt*(qp(k-1)-qb(k-1))/...
(const_mud_pipe2+const_cem_pipe2+const_spacer_pipe2));
Pdha(k) = max(0,Pdha(k-1) + dt*(B_mud/V_ta)*(qb(k-1) - qout(k-1)));
Ppump(k) = max(0,Pp(k)+Fsurf(k));
elseif statec(k) == 3
Vcp_c(k) = Vcp_c(k-1);
Vsp_c1(k) = Vsp_c1(k-1);
Vmp_c1(k) = max(0,Vmp1(k-1) - dt*(qb(k-1)));
Vsp_c2(k) = Vsp2(k-1) + dt*(qp(k-1));
Vma(k) = V_ta;
Vma_c(k) = V_ta;
const_mud_pipe1 = Vmp_c1(k)/B_mud;
const_spacer1_pipe1 = Vsp_c1(k)/B_spacer;
const_spacer2_pipe1 = Vsp_c2(k)/B_spacer2;
const_cem_pipe1 = Vcp_c(k)/B_cem;
deltaP(k) = -(V_tp-hd(k-1)*Ap-Vmp_c1(k)-Vsp_c1(k)-Vsp_c2(k)-Vcp_c(k))/...
(const_mud_pipe1+const_spacer1_pipe1+ const_spacer2_pipe1+...
const_cem_pipe1);
Vmp1(k) = Vmp_c1(k)*(1-deltaP(k)/B_mud);
Vsp2 (k) = Vsp_c2(k)*(1-deltaP(k)/B_spacer2);
Vsp1(k) = Vsp_c1(k)*(1-deltaP(k)/B_spacer);
Vcp(k) = Vcp_c(k)*(1-deltaP(k)/B_cem);
Lcp(k) = (Vcp(k)/Ap);
Lcpf(k) = (Vcp(k)/Ap)+(Vsp2(k)+Vmp2(k))/Ap+hd(k-1);
Lsp1(k) = Vsp1(k)/Ap;
Lsp2(k) = Vsp2(k)/Ap;
Lspf(k) = Lcpf(k) + Lsp1(k);
Lspf2(k) = (Vsp2(k)/Ap)+hd(k-1)+(Vmp2(k)/Ap);
const_mud_pipe2 = Vmp1(k)/B_mud;
const_mud2_pipe2 = Vmp2(k)/B_mud;
const_spacer1_pipe2 = Vsp1(k)/B_spacer;
const_spacer2_pipe2 = Vsp2(k)/B_spacer2;
const_cem_pipe2 = Vcp(k)/B_cem;
Fsurf(k) = Ffactor*unified_friction_pipe(Surfd_i,qp(k-1),L_surf,PV_spacer2,YP_SI_spacer2,tau_y_spacer2,rho_spacer2);
Pp(k) = max(0,Pp(k-1) + dt*(qp(k-1)-qb(k-1))/(const_mud_pipe2+ ...
const_mud2_pipe2+...
const_spacer1_pipe2+const_spacer2_pipe2+const_cem_pipe2));
Pdha(k) = Pdha(k-1) + dt*(B_mud/V_ta)*(qb(k-1)-qout(k-1));
Ppump(k) = max(0,Pp(k)+Fsurf(k));
elseif statec(k) == 4
Vcp_c(k) = Vcp_c(k-1);
Vsp_c1(k) = Vsp_c1(k-1);
Vsp_c2(k) = Vsp_c2(k-1);
Vmp_c1(k) = max(0,Vmp1(k-1) - dt*(qb(k-1)));
Vmp_c2(k) = Vmp2(k-1) + dt*(qp(k-1));
Vma(k) = V_ta;
Vma_c(k) = V_ta;
const_mud_pipe1 = Vmp_c1(k)/B_mud;
const_mud2_pipe1 = Vmp_c2(k)/B_mud;
const_spacer1_pipe1 = Vsp_c1(k)/B_spacer;
const_spacer2_pipe1 = Vsp_c2(k)/B_spacer2;
const_cem_pipe1 = Vcp_c(k)/B_cem;
deltaP(k) = -(V_tp-hd(k-1)*Ap-Vmp_c1(k)-Vsp_c1(k)-Vsp_c2(k)-Vcp_c(k)-...
Vmp_c2(k))/...
(const_mud_pipe1+const_mud2_pipe1+const_spacer1_pipe1+...
const_spacer2_pipe1+ const_cem_pipe1);
Vmp1(k) = Vmp_c1(k)*(1-deltaP(k)/B_mud);
Vsp2 (k) = Vsp_c2(k)*(1-deltaP(k)/B_spacer2);
Vsp1(k) = Vsp_c1(k)*(1-deltaP(k)/B_spacer);
Vmp2(k) = Vmp_c2(k)*(1-deltaP(k)/B_mud);
Vcp(k) = Vcp_c(k)*(1-deltaP(k)/B_cem);
Lcp(k) = (Vcp(k)/Ap);
Lcpf(k) = (Vcp(k)/Ap)+(Vsp2(k)/Ap)+hd(k-1)+(Vmp2(k)/Ap);
Lsp1(k) = Vsp1(k)/Ap;
Lsp2(k) = Vsp2(k)/Ap;
Lspf(k) = Lcpf(k) + Lsp1(k);
Lspf2 (k) =(Vsp2(k)/Ap)+hd(k-1)+(Vmp2(k)/Ap); %+(Vmp2(k)/Ap)
const_mud_pipe2 = Vmp1(k)/B_mud;
const_spacer1_pipe2 = Vsp1(k)/B_spacer;
const_spacer2_pipe2 = Vsp2(k)/B_spacer2;
const_cem_pipe2 = Vcp(k)/B_cem;
Fsurf(k) = Ffactor*unified_friction_pipe(Surfd_i,qp(k-1),L_surf,PV_mud,YP_SI_mud,tau_y_mud,rho_mud);
Pp(k) = max(0,Pp(k-1) + dt*(qp(k-1)-qb(k-1))/(const_mud_pipe2+ ...
const_spacer1_pipe2+const_spacer2_pipe2+const_cem_pipe2));
Pdha(k) = Pdha(k-1) + dt*(B_mud/V_ta)*(qb(k-1)-qout(k-1));
Ppump(k) = max(0,Pp(k)+Fsurf(k));
elseif statec(k) == 5
% Pipe
Vcp_c(k) = Vcp_c(k-1);
Vsp_c2(k) = Vsp_c2(k-1);
Vmp_c2(k) = Vmp2(k-1) + dt*(qp(k-1));
Vsp_c1(k) = max(0,Vsp1(k-1) - dt*(qb(k-1)));
const_mud_pipe1 = Vmp_c2(k)/B_mud;
const_spacer1_pipe1 = Vsp_c1(k)/B_spacer;
const_spacer2_pipe1 = Vsp_c2(k)/B_spacer2;
const_cem_pipe1 = Vcp_c(k)/B_cem;
deltaP(k) = -(V_tp-hd(k-1)*Ap-Vmp_c2(k)-Vcp_c(k)-Vsp_c1(k)-Vsp_c2(k))/...
(const_mud_pipe1+const_spacer1_pipe1+const_spacer2_pipe1+...
const_cem_pipe1);
Vmp1(k) = 0;
Vmp_c1(k) = 0;
(const_mud_pipe1+const_spacer1_pipe1+const_spacer2_pipe1+const_cem_pipe1);
Vmp2(k) = Vmp_c2(k)*(1-deltaP(k)/B_mud);
Vcp(k) = max(0,Vcp_c(k)*(1-deltaP(k)/B_cem));
Vsp1(k) = max(0,Vsp_c1(k)*(1-deltaP(k)/B_spacer));
Vsp2(k) = max(0,Vsp_c2(k)*(1-deltaP(k)/B_spacer2));
Lcp(k) = (Vcp(k)/Ap);
Lcpf(k) = (Vcp(k)/Ap)+(Vsp2(k)/Ap)+hd(k-1)+(Vmp2(k)/Ap);
Lsp1(k) = Vsp1(k)/Ap;
Lsp2(k) = Vsp2(k)/Ap;
Lspf(k) = Lcpf(k) + Lsp1(k);
Lspf2(k) =(Vsp2(k)/Ap)+hd(k-1)+(Vmp2(k)/Ap);
const_mud_pipe2 = Vmp2(k)/B_mud;
const_spacer1_pipe2 = Vsp1(k)/B_spacer;
const_spacer2_pipe2 = Vsp2(k)/B_spacer2;
const_cem_pipe2 = Vcp(k)/B_cem;
Fsurf(k) = Ffactor*unified_friction_pipe(Surfd_i,qp(k-1),L_surf,PV_mud,YP_SI_mud,tau_y_mud,rho_mud);
Pp(k) = max(0,Pp(k-1) + dt*(qp(k-1)-qb(k-1))/(const_mud_pipe2+...
const_cem_pipe2+const_spacer2_pipe2+const_spacer1_pipe2));
Ppump(k) = max(0,Pp(k)+Fsurf(k));
% Annulus
Vma_c(k) = Vma(k-1) - dt*(qout(k-1));
Vsa_c(k) = Vsa(k-1) + dt*(qb(k-1));
const_spacer_an1 = Vsa_c(k)/B_spacer;
const_mud_an1 = Vma_c(k)/B_mud;
deltaP(k) = -(V_ta-Vma_c(k)-Vsa_c(k))/(const_mud_an1+const_spacer_an1);
Vma(k) = Vma_c(k)*(1-deltaP(k)/B_mud);
Vsa(k) = Vsa_c(k)*(1-deltaP(k)/B_spacer);
Lsa(k) = Vsa(k)/Aa;
const_mud_an2 = Vma(k)/B_mud;
const_spacer_an2 = Vsa(k)/B_spacer;
Pdha(k) = Pdha(k-1) + dt*(qb(k-1)-qout(k-1))/(const_mud_an2+const_spacer_an2);
elseif statec(k) == 6
% Pipe
Vcp_c(k) = max(0,Vcp(k-1)-dt*(qb(k-1)));
Vsp_c2(k) = Vsp_c2(k-1);
Vmp_c2(k) = Vmp2(k-1) + dt*(qp(k-1));
const_mud_pipe1 = Vmp_c2(k)/B_mud;
const_spacer2_pipe1 = Vsp_c2(k)/B_spacer2;
const_cem_pipe1 = Vcp_c(k)/B_cem;
deltaP(k) = -(V_tp-hd(k-1)*Ap-Vmp_c2(k)-Vcp_c(k)-Vsp_c2(k))/...
( const_mud_pipe1+const_spacer2_pipe1+const_cem_pipe1);
Vmp1(k) = 0;
Vmp_c1(k) = 0;
Vsp_c1(k) = 0;
Vsp1(k) = 0;
Vmp2(k) = Vmp_c2(k)*(1-deltaP(k)/B_mud);
Vcp(k) = max(0,Vcp_c(k)*(1-deltaP(k)/B_cem));
Vsp2(k) = max(0,Vsp_c2(k)*(1-deltaP(k)/B_spacer2));
Lcp(k) = (Vcp(k)/Ap);
Lsp1(k) = 0;
Lsp2(k) = Vsp2(k)/Ap;
Lspf2(k) =(Vsp2(k)/Ap)+hd(k-1)+(Vmp2(k)/Ap);
const_mud_pipe2 = Vmp2(k)/B_mud;
const_spacer2_pipe2 = Vsp2(k)/B_spacer2;
const_cem_pipe2 = Vcp(k)/B_cem;
Fsurf(k) = Ffactor*unified_friction_pipe(Surfd_i,qp(k-1),L_surf,PV_mud,YP_SI_mud,tau_y_mud,rho_mud);
Pp(k) = max(0,Pp(k-1) + dt*(qp(k-1)-qb(k-1))/(const_mud_pipe2+...
const_cem_pipe2+const_spacer2_pipe2));
Ppump(k) = max(0,Pp(k)+Fsurf(k));
% Annulus
Vma_c(k) = Vma(k-1) - dt*(qout(k-1));
Vsa_c(k) = Vsa_c(k-1);
Vca_c(k) = Vca(k-1) + dt*(qb(k-1));
const_spacer_an1 = Vsa_c(k)/B_spacer;
const_mud_an1 = Vma_c(k)/B_mud;
const_cem_an1 = Vca_c(k)/B_cem;
deltaP(k) = -(V_ta-Vma_c(k)-Vsa_c(k)-...
Vca_c(k))/(const_mud_an1+const_spacer_an1+const_cem_an1);
Vma(k) = Vma_c(k)*(1-deltaP(k)/B_mud);
Vsa(k) = Vsa_c(k)*(1-deltaP(k)/B_spacer);
Vca(k) = Vca_c(k)*(1-deltaP(k)/B_cem);
Lsa(k) = Vsa(k)/Aa;
Lca(k) = Vca(k)/Aa;
const_mud_an2 = Vma(k)/B_mud;
const_spacer_an2 = Vsa(k)/B_spacer;
const_cem_an2 = Vca(k)/B_cem;
Pdha(k) = Pdha(k-1) + dt*(qb(k-1)-qout(k-...
1))/(const_mud_an2+const_spacer_an2+const_cem_an2);
elseif statec(k) == 7 %Pipe
Vsp_c2(k) = max(0,Vsp2(k-1)-dt*(qb(k-1)));
Vmp_c2(k) = Vmp2(k-1) + dt*(qp(k-1));
const_mud_pipe1 = Vmp_c2(k)/B_mud;
const_spacer2_pipe1 = Vsp_c2(k)/B_spacer2;
deltaP(k) = -(V_tp-hd(k-1)*Ap-Vmp_c2(k)-Vsp_c2(k))/( const_mud_pipe1+const_spacer2_pipe1);
Vmp1(k) = 0;
Vmp_c1(k) = 0;
Vsp_c1(k) = 0;
Vsp1(k) = 0;
Vcp (k) = 0 ;
Vcp_c(k) = 0 ;
Vmp2(k) = Vmp_c2(k)*(1-deltaP(k)/B_mud);
Vsp2(k) = max(0,Vsp_c2(k)*(1-deltaP(k)/B_spacer));
const_mud_pipe2 = Vmp2(k)/B_mud;
const_spacer2_pipe2 = Vsp2(k)/B_spacer2;
Lcp(k) = 0;
Lsp1(k) = 0;
Lsp2(k) = Vsp2(k)/Ap;
Fsurf(k) = unified_friction_pipe(Surfd_i,qp(k-1),L_surf,PV_mud,YP_SI_mud,tau_y_mud,rho_mud);
Pp(k) = max(0,Pp(k-1) + dt*(qp(k-1)-qb(k-1))/(const_mud_pipe2+const_spacer2_pipe2));
Ppump(k) = max(0,Pp(k)+Fsurf(k));
%Annulus
Vma_c(k) = Vma(k-1) - dt*(qout(k-1));
Vsa_c(k) = Vsa_c(k-1);
Vca_c(k) = Vca_c(k-1);
Vsa_c2(k) = Vsa2(k-1)+ dt*(qb(k-1));
const_spacer_an1 = Vsa_c(k)/B_spacer;
const_mud_an1 = Vma_c(k)/B_mud;
const_cem_an1 = Vca_c(k)/B_cem;
const_spacer2_an1= Vsa_c2 (k)/B_spacer2;
deltaP(k) = -(V_ta-Vma_c(k)-Vsa_c(k)-Vca_c(k)-Vsa_c2(k))/(const_mud_an1+const_spacer_an1+const_cem_an1+const_spacer2_an1);
Vma(k) = Vma_c(k)*(1-deltaP(k)/B_mud);
Vsa(k) = Vsa_c(k)*(1-deltaP(k)/B_spacer);
Vca(k) = Vca_c(k)*(1-deltaP(k)/B_cem);
Vsa2(k) = Vsa_c2(k)*(1-deltaP(k)/B_spacer2);
Lsa(k) = Vsa(k)/Aa;
Lca(k) = Vca(k)/Aa;
Lsa2(k) = Vsa2(k)/Aa;
const_mud_an2 = Vma(k)/B_mud;
const_spacer_an2 = Vsa(k)/B_spacer;
const_cem_an2 = Vca(k)/B_cem;
const_spacer2_an2 = Vsa2(k)/B_spacer2;
Pdha(k) = Pdha(k-1) + dt*(qb(k-1)-qout(k-1))/(const_mud_an2+const_spacer_an2+const_cem_an2+const_spacer2_an2);
else
% qp(k) = max(0, qp(k-1) - dt*(qp_mud/tramp));
%Pipe
Vsp_c2(k) = max(0,Vsp2(k-1)-dt*(qb(k-1)));
Vmp_c2(k) = Vmp2(k-1) + dt*(qp(k-1));
const_mud_pipe1 = Vmp_c2(k)/B_mud;
const_spacer2_pipe1 = Vsp_c2(k)/B_spacer2;
deltaP(k) = -(V_tp-hd(k-1)*Ap-Vmp_c2(k)-Vsp_c2(k))/( const_mud_pipe1+const_spacer2_pipe1);
Vmp1(k) = 0;
Vmp_c1(k) = 0;
Vsp_c1(k) = 0;
Vsp1(k) = 0;
Vcp (k) = 0 ;
Vcp_c(k) = 0 ;
Vmp2(k) = Vmp_c2(k)*(1-deltaP(k)/B_mud);
Vsp2(k) = max(0,Vsp_c2(k)*(1-deltaP(k)/B_spacer));
const_mud_pipe2 = Vmp2(k)/B_mud;
const_spacer2_pipe2 = Vsp2(k)/B_spacer2;
Lcp(k) = 0;
Lsp1(k) = 0;
Lsp2(k) = Vsp2(k)/Ap;
Fsurf(k) = unified_friction_pipe(Surfd_i,qp(k-1),L_surf,PV_mud,YP_SI_mud,tau_y_mud,rho_mud);
Pp(k) = max(0,Pp(k-1) + dt*(qp(k-1)-qb(k-1))/(const_mud_pipe2+const_spacer2_pipe2));
Ppump(k) = max(0,Pp(k)+Fsurf(k));
%Annulus
Vma_c(k) = Vma(k-1) - dt*(qout(k-1));
Vsa_c(k) = Vsa_c(k-1);
Vca_c(k) = Vca_c(k-1);
Vsa_c2(k) = Vsa2(k-1)+ dt*(qb(k-1));
const_spacer_an1 = Vsa_c(k)/B_spacer;
const_mud_an1 = Vma_c(k)/B_mud;
const_cem_an1 = Vca_c(k)/B_cem;
const_spacer2_an1= Vsa_c2 (k)/B_spacer2;
deltaP(k) = -(V_ta-Vma_c(k)-Vsa_c(k)-Vca_c(k)-Vsa_c2(k))/(const_mud_an1+const_spacer_an1+const_cem_an1+const_spacer2_an1);
Vma(k) = Vma_c(k)*(1-deltaP(k)/B_mud);
Vsa(k) = Vsa_c(k)*(1-deltaP(k)/B_spacer);
Vca(k) = Vca_c(k)*(1-deltaP(k)/B_cem);
Vsa2(k) = Vsa_c2(k)*(1-deltaP(k)/B_spacer2);
Lsa(k) = Vsa(k)/Aa;
Lca(k) = Vca(k)/Aa;
Lsa2(k) = Vsa2(k)/Aa;
const_mud_an2 = Vma(k)/B_mud;
const_spacer_an2 = Vsa(k)/B_spacer;
const_cem_an2 = Vca(k)/B_cem;
const_spacer2_an2 = Vsa2(k)/B_spacer2;
Pdha(k) = Pdha(k-1) + dt*(qb(k-1)-qout(k-1))/(const_mud_an2+const_spacer_an2+const_cem_an2+const_spacer2_an2);
end
if hd(k-1) ~= 0
Pp(k) = 0;
end
%% Dynamics for the total wellbore
qb(k) = max(0,qb(k-1) + dt/(Mp(k-1)+Ma(k-1))*((Pp(k-1) - P0 - Fp(k-1)-Fa(k-1) +...
Gp(k-1)-Ga(k-1))));
qavg_pipe(k) = max(0,qavg_pipe(k-1) + dt/Mp(k-1)*(Pp(k-1) - Fp(k-1) - Pdha(k-1) + Gp(k-1)));
qout(k) = max(0,qout(k-1) + dt/(Ma(k-1))*(Pdha(k-1) - Fa(k-1) - P0 - Ga(k-1)));
if Pp(k) ==0 && ((Gp(k-1)-Ga(k-1)-Fa(k-1)-Fp(k-1))) ~= 0
hd(k) = max(0,hd(k-1) + (dt/Ap)*(qb(k-1)-qp(k-1)));
zz(k)=5; % to know if this condition is work or not
else
hd(k) = hd(k-1);
zz(k)=-5;
end
% if k >= 97*6000 && k<= 122*6000
% Pp(k) = 165/14.5*10^5 ;
% end
%
%% Update variables
Lsp(k) = Lsp1(k) + Lsp2(k);
Lma(k) = L_SI - Lca(k) - Lsa(k) - Lsa2(k);
Lmp(k) = L_SI - Lcp(k) - hd(k)-Lsp(k);
q_avgp = qavg_pipe(k);
q_avga = qout(k);
Mp(k) = (rho_mud*Lmp(k)/(Ap) + rho_cem*Lcp(k)/(Ap) + rho_spacer*Lsp1(k)/(Ap)+ rho_spacer2*Lsp2(k)/(Ap));
Ma(k) = (rho_mud*Lma(k)/(Aa) + rho_cem*Lca(k)/(Aa) + rho_spacer*Lsa(k)/(Aa) + rho_spacer2*Lsa2(k)/(Ap));
Gp(k) = rho_mud*g*Lmp(k) + rho_cem*g*Lcp(k) + rho_spacer*g*Lsp1(k) + rho_spacer2*g*Lsp2(k);
Ga(k) = rho_mud*g*Lma(k) + rho_cem*g*Lca(k) + rho_spacer*g*Lsa(k) + rho_spacer2*g*Lsa2(k);
Fp(k) = Ffactor*unified_friction_pipe(Sd_i,q_avgp,Lmp(k),PV_mud,YP_SI_mud,tau_y_mud,rho_mud) +...
unified_friction_pipe(Sd_i,q_avgp,Lsp1(k),PV_spacer,YP_SI_spacer,tau_y_spacer,rho_spacer)+...
unified_friction_pipe(Sd_i,q_avgp,Lcp(k),PV_cem,YP_SI_cem,tau_y_cem,rho_cem)+...
unified_friction_pipe(Sd_i,q_avgp,Lsp2(k),PV_spacer2,YP_SI_spacer2,tau_y_spacer2,rho_spacer2);
Fa(k) = Ffactor*unified_friction_ann(Cd_o,Od,q_avga,Lma(k),PV_mud,YP_SI_mud,tau_y_mud,rho_mud) +...
unified_friction_ann(Cd_o,Od,q_avga,Lsa(k),PV_spacer,YP_SI_spacer,tau_y_spacer,rho_spacer)+...
unified_friction_ann(Cd_o,Od,q_avga,Lca(k),PV_cem,YP_SI_cem,tau_y_cem,rho_cem)+...
unified_friction_ann(Cd_o,Od,q_avga,Lsa2(k),PV_spacer2,YP_SI_spacer2,tau_y_spacer2,rho_spacer2);
if qp(k) == 0
Ppump(k) = 0;
end
end
% %% The simulation is finished. The simulation results can then be plotted as desired.
% % %plot realpump pressure data %%
TT=readtable('data13_2.xlsx'); %excel filename
tpr=TT(:,1);
ppr=TT(:,5); %139*1
%% Same length
Px=Ppump(1:2400:N); %139 cases
Pxt=transpose(Px); %138*1xl
%ecd=Pdha/14.5*10.2/L_SI;
ecd= rho_mud +(Pdha*10.2/L_SI);
%plotres
tm=t/60;
figure(1)
plot(tm,qp/lpm2m3s)
grid
xlabel('Time [min]')
ylabel('Flow [lpm]')
%legend('Pump','Bit','Out')
figure(2)
%plot(tm,Ppump*1e-6*14.5,tm,Pp*1e-6*14.5)
plot(tm,Ppump*1e-6*14.5)
grid
xlabel('Time [min]')
ylabel('Pressure [psi]')
%legend('Pump','Top of pipe')
legend('Pump pressure')
figure(3)
plot(tm,hd)
grid
xlabel('Time [min]')
ylabel('Free fall [m]')
figure(4)
plot(tm,Fsurf*1e-6*14.5,tm,Fp*1e-6*14.5,tm,Fa*1e-6*14.5)
grid
xlabel('Time [min]')
ylabel('Friction pressure [psi]')
legend('Surface','Pipe','Annulus')
figure(5)
plot(tm,ecd)
grid
xlabel('Time [min]')
ylabel('ECD [kg/l]')
  1 Comment
Torsten
Torsten on 27 May 2023
Edited: Torsten on 27 May 2023
How should we be able to correct this code ?
Excel data and functions are missing to execute it.
And if results are not as expected, we can help less than ever.

Sign in to comment.

Sign in to answer this question.

Answers (2)

dpb
dpb on 27 May 2023
Moved: dpb on 27 May 2023
It's so well documented, should be a breeze...
The main problem appears to be the missing function giving the error
Unrecognized function or variable 'unified_friction_pipe'.
If as it appears you didn't write it, return to the place from whence you got this part and download the rest.
  2 Comments
noobie
noobie on 5 Jun 2023
I have the function it works fine and the code does not propt errors jut giving the the wrong results
Walter Roberson
Walter Roberson on 5 Jun 2023
I would not expect any of the volunteers to go through the effort of digging the function out of the book I listed, and run your code, and immediately know not only what results you see on your system but also know what the "right" results are and debug the code for you

Sign in to comment.

Walter Roberson
Walter Roberson on 27 May 2023
The missing code for unified_friction_pipe is on page 100 of
Managed Pressure Cementing
Simulations of Pressure and Flow Dynamics During Cementing Using Applied Back- Pressure and Dual Gradient
by
Øystein Seglem Bekken
Erik Havnen Ullsfoss

Categories

Find more on Functions in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!