%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Ane Cecilie Duus
%Master thesis 06/06/13

%Script for calculating and plotting:
%   1. Bottom hole pressure 
%   2. Pressure profile in Pipe
%   3. Density profile in Pipe
%   4. Pressure drop across valve

%%


load press_drop_valve;      %Pressure drop profile (valve) from OLGA
load bottom_press;          %Bottom pressure from OLGA

figure
hold on


P_res=205;      %   Reservoir pressure [bara]
T_res= 63;      %   Reservoir temperature [C]

P_bh=P_res;     %   Initial value for bottom hole pressure
k=0;            %   Define start value for k=0

for i= 0:0.01:1     %   For valve opening (vo) = 0:1

    
    vo=i;           %   Define i = vo (valve opening)
    k=k+1;
    
    
    
 %%%%% FLOW ESTIMATES FROM OLGA  %%%%%
 
 
 % Flow-profiles fitted to data from OLGA-simulations. 
 
q_gas1= (3.05*10^(4))*(1-exp(-40*vo));
Q_GAS(k)=q_gas1;                    
 
q_oil=(155)*(1-exp(-40*vo));
Q_OIL(k)=q_oil;

q_water=(287)*(1-exp(-40*vo));
Q_WATER(k)=q_water;

q_gasm=q_gas1*(1/P_bh)

q_tot=(q_gasm+q_oil+q_water)/60/60/24;      %   Total flow [m^3/s]




%%%% CALCULATION OF BOTTOMHOLE PRESSURE %%%%
mu=1050;                 %  Viscosity
kh=0.0999;               %  Permability(k) x formation thickness(h)

   
       %%%%%%%%%%% DARCY'S LAW %%%%%%%%%%%

P_bh= 250-(q_tot*mu*(log(1000/0.1))/(2*3.14*kh))    % Bottom hole pressure [bar]



PW(k)=P_bh;              %Vector for bottom hole pressures for vo=0.01-1




%%%% CALCUALTION OF PRESSURE & DENSITY PROFILE IN PIPE %%%%

par.R = 8.314;      % Gas Constant
par.T = 336.15;     % Reservoir temperature [K]
par.Mg = 0.0195;    % Molar weight of natural gas [kg/mol]
par.rho_o = 945;    % Density of oil [kg/m^3]
par.g  = 9.81;
par.D = 0.12;       % m


%%% VOLUME ESTIMATIONS %%%

% Volumes estimated from data obtained by simulations in OLGA
par.vtot_0= -3119.2*(vo^6) + 11356*(vo^5) - 16537*(vo^4) + 12292*(vo^3) - 4918.2*(vo^2) + 1013.3*(vo) + 566.25;
par.vl= -1902.3*(vo^6) + 6923.1*(vo^5) - 10077*(vo^4) + 7485.1*(vo^3) - 2992.1*(vo^2) + 615.57*(vo) + 425.01;


%%% INITIAL PRESSURE IN BOTTOM OF PIPELINE %%%

par.P_0= P_bh*10e4;   % [Pa]


%%% CALCULATION OF GAS VOLUME %%%

vg_0= par.vtot_0-par.vl;
par.n= par.P_0*vg_0/(par.R*par.T*0.85);

p_0= par.P_0;


%%% CALCULATION OF PRESSURE PROFILE IN THE PIPE
% Ode23s calculates pressureprofile from bottom hole (inflow, at 1300-26
% meters depth) to valve at 1300 - 1300 meters depth.

p = ode23s(@(z,p) pressure_drop_model(z,p,par), [26 1300],p_0);

P=p.y';

%%% COMPRESSIBILITY FACTOR FROM BOTTOM TO TOP %%%

Z = 0.85+(0.15)/1300*p.x';      % Linear increas of Z-factor.



%%% DENSITY CALCULATIONS %%%

vg=par.n*par.R*par.T./P;                             %  Calculation og gas-volume
vl=par.vl*ones(size(vg));                            %  Calculation of liquid-volume
vtot =vg+vl;                                         %  w_g/par.w;
  
 rho_g =P.*par.Mg./(par.R.*par.T.*Z);                %  Gas-density using compressibilty-factor   
 
 rho_o=par.rho_o*ones(size(vtot));                   %  Density profile   
 rho = (vg./vtot).*rho_g + (1-(vg./vtot)).*rho_o;    %  Calculation of mixed density   [kg/m^3]  

 
 % Interpolation to obtain initial wellhead pressure
 
x0 = 15:2:1300;
y0=INTERP1(p.x,p.y,x0);         % Interpoletes pressure 
P_wh=y0(559)*10e-6              % Pressure before valve 
 
P_WH(k)=P_wh;


%%%% CALCULTION OF PRESSUREDROP OVER VALVE    %%%%

A1=0.0046; % Orfice cross-sectional area [m^2]
Q=0.0036; % Total mass flow [m3/s]
Cd=0.84; % Discharge coefficient



% Interpolation to obtain density at well head

x0 = 15:2:1300;
rho0=INTERP1(p.x,rho,x0);   % Interpoletes pressure 
rho_wh=rho0(559);           % Density before valve 



%%% CALCULATIONS OF PRESSURE DROP ACROSS VALVE %%%

x=vo;           
fz = 0.13*(1-exp(-10*x));   % Charecteristic Valve Equation 
FZ(k)=fz;
 


 %Bernoulli
 
 %dP_valve= 0.5*rho_wh*((q_tot/A1)^2)*((1/(vo*Cd))^2-1)*10e-6 ; 
 
 
 %Bernoulli+dp_min
 
 %dP_valve= 0.5*rho_wh*((q_tot/A1)^2)*((1/(vo*Cd))^2 + (1/vo)^2 )*10e-6 + 1.3;  
 
 
 %Valve Equation
 
 dP_valve= ((q_tot/((fz)*A1*Cd))^2)*rho_wh*10e-6 ;
 
 
 
 P_out=P_wh-dP_valve;       % Pressure at valve outlet
  
 dp_VALVE(k)=dP_valve;      % Pressure drop valve
 P_OUT(k)=P_out;

 
end



%%%%    PLOTS    %%%%%



%%% Inflow vs. vo Plots %%%

%load inflow_vo

%Gas
%plot(0:0.01:1,Q_GAS,'r','linewidth',2)
%hold on
% plot(inflow_vo(:,1),inflow_vo(:,2),'linewidth',2)
% ylabel 'Gas Flow [bar]'
% title 'Inflow Gas'

%Oil
% plot(0:0.01:1,Q_OIL,'r','linewidth',2)
% hold on
% plot(inflow_vo(:,1),inflow_vo(:,3),'linewidth',2)
% ylabel 'Oil Flow [bar]'
% title 'Inflow Oil'

%Water
% plot(0:0.01:1,Q_WATER,'r','linewidth',2)
% hold on
% plot(inflow_vo(:,1),inflow_vo(:,4),'linewidth',2)
% legend('Estimation','OLGA')
% ylabel 'Gas Flow [bar]'
% title 'Inflow Gas'
% xlabel 'Valve Opening'

%PW=PW(1,2:102);

%plot(0:0.01:1,dp_VALVE,'r','linewidth',2)
%plot(0:0.01:1,dp_VALVE,'r','linewidth',2)
%plot(0:0.01:1,PW,'r','linewidth',2)
%plot(0:0.01:1,dens,'r','linewidth',2)
%hold on
plot(bottom_press(:,1),bottom_press(:,2),'linewidth',2)
%plot(press_drop_valve(:,1),press_drop_valve(:,2),'linewidth',2)
%legend('Estimation','OLGA')
xlabel 'Valve Opening'
ylabel 'Pressure [bar]'
title 'Bottom Hole Pressure OLGA'
axis([0 1 120 250])
%title 'Bottom Pressure'
%title 'Wellhead Pressure'
%title 'Pressure after Valve(Cd=50)'
%title 'Pressure Drop Across Valve (Valve eq)'






%%% Flow Plots %%%

% load inflow
% plot(inflow(:,1),inflow(:,2)/1000,'r','linewidth',2)
% hold on
% plot(inflow(:,1),inflow(:,3),'b','linewidth',2)
% hold on
% plot(inflow(:,1),inflow(:,4),'g','linewidth',2)
% hold on
% legend('Gas Flow','Oil Flow','Water Flow')
% xlabel 'Pressure'
% ylabel 'Flow []'
% title 'Inflow From Reservoir'
% hold off

%%% q plot %%%

%  load fz
%  plot(fz(:,1),fz(:,2),'b','linewidth',2)
%  xlabel 'Valve Opening'
%  ylabel 'Flow [% of max flow]'
%  title 'Flow Characteristics for Valve'


%%% fz plot %%%

% load fz_curve2
% plot(0:0.01:1,FZ,'r','linewidth',2)
% hold on
%  plot(fz_curve2(:,1),fz_curve2(:,7),'b','linewidth',2)
%  legend('Estimation','OLGA')
%  xlabel 'Valve Opening'
%  ylabel 'f(z)'
%  axis([0 1 0 0.3])
%  title 'Characteristic Function f(z) for Valve OLGA'
 

% hold all
% 
%plot(0:0.01:1,P_WH)
% hold all
% 
%plot(0:0.01:1,P_OUT)
% %axis([20 105 0 2000])

hold off

 %P_BH((i*100)+1,1)=P_bh;
%P_WH((i*100)+1,1)=P_wh;
%P_AW((i*100)+1,1)=P_out;