% Simplified model of twophase pipeline with slug behavior
% Gisle Otto Eikrem
% 15.05.2004

m_gas_riser = m(1,i); 
m_oil_riser = m(2,i);
m_gas_pipeline = m(3,i); 
m_oil_pipeline = m(4,i);

rho_water = 1000; %kg/m3, water
  
%------------------------------------------------------------------------  

% Calculation of pressure in pipeline section
  
V_oil_pipeline = m_oil_pipeline/rho_oil; %Assumes known density of the oil
  
V_gas_pipeline = (V_total_pipeline - V_oil_pipeline); %The volume of gas in the pipeline is given when the volume of oil is known.
  
p_top_pipeline_bara(i) = m_gas_pipeline*R1*T_pipeline/(V_gas_pipeline*M); %Calculate the pressure at the inlet of the pipeline

rho_gas_top_pipeline(i) = m_gas_pipeline / V_gas_pipeline; % density of gas in pipeline

rho_m_pipeline(i) = (m_gas_pipeline + m_oil_pipeline)/(pipeline_length*A_pipeline); % Mixture density in pipeline 
  
delta_p_pipeline(i) = rho_m_pipeline(i)*g*pipeline_length*sin(alpha_pipeline)*conv_to_bara; % bara, pressure drop in tubing down to the bottom of the well   

delta_p_gas(i) = rho_gas_top_pipeline(i)*g*(V_gas_pipeline/A_pipeline)*sin(alpha_pipeline)*conv_to_bara;

p_bottom_pipeline_bara_oil(i) = p_top_pipeline_bara(i) + delta_p_pipeline(i); %bara, Pressure in bottom of pipeline

p_bottom_pipeline_bara_gas(i) = p_top_pipeline_bara(i) + delta_p_gas(i); %bara, Pressure in bottom of pipeline

p_gas_flow_noflow(i) = p_top_pipeline_bara(i) + rho_m_pipeline(i) * g * pipeline_length * sin(alpha_pipeline);
  
% Calculation of pressure in riser section
  
V_oil_riser = m_oil_riser/rho_oil; %Assumes known density of the oil
  
V_gas_riser = (V_total_riser - V_oil_riser); %The volume of gas in the riser is given when the volume of oil is known.

m_gas_riser_min(i) = p_sep_bara * V_gas_riser * M / (R1 * T_riser); % Calculation of the lower limit of gas in the riser, maintained by backflow from the separator
  
p_top_riser_bara(i) = m_gas_riser*R1*T_riser/(V_gas_riser*M); %Calculate the pressure at the inlet of the riser
  
rho_m_riser(i) = (m_gas_riser + m_oil_riser)/(riser_length*A_riser); % Mixture density in riser 
  
delta_p_riser(i) = rho_m_riser(i)*g*riser_length*sin(alpha_riser)*conv_to_bara; % bara, pressure drop in tubing down to the bottom of the well   
  
p_bottom_riser_bara(i) = p_top_riser_bara(i) + delta_p_riser(i); %bara, Pressure in bottom of riser

p_bottom_riser_onlyoil_bara(i) = p_top_riser_bara(i) + (rho_oil * g * riser_length * sin(alpha_riser))*conv_to_bara;


% dynamic GOR inflow or constant
if strcmp(GOR_type,'constant')
        X_gas(i)=X_gas_nom;
        X_oil(i)=1-X_gas(i);
elseif strcmp(GOR_type,'linear')
       X_gas(i)=X_gas_nom + X_gas_nom*0.1*(p_top_pipeline_bara(i) - 210);   
        X_oil(i)=1-X_gas(i);
end

% Inflow to the pipeline

if flow_type == 'constant '
    
w_oil_in(i) = inflow_oil;
w_gas_in(i) = inflow_gas;

elseif flow_type == 'linear   ';

w_oil_in(i) = X_oil(i)*1.0107*1.0245*prodi*(p_reservoir - p_top_pipeline_bara(i));
w_gas_in(i) = (X_gas(i)*w_oil_in(i))./(1-X_gas(i));    
    
elseif flow_type == 'reservoir';

reservoir_model;
w_oil_in(i)=X_oil(i)*w_tot_in(i);
w_gas_in(i) = X_gas(i)*w_tot_in(i); %(X_gas(i)*w_oil_in(i))./(1-X_gas(i));    
end



% Flow across bottom point
G_pipeline_bottom_oil(i) = rho_oil/rho_water;  
G_pipeline_bottom_gas(i) = 11/rho_water;

q_bottom_oil(i) = (p_bottom_pipeline_bara_oil(i) > p_bottom_riser_bara(i))*K_r_oil*sqrt((p_bottom_pipeline_bara_oil(i) - p_bottom_riser_bara(i))/G_pipeline_bottom_oil(i)); %Flow of oil from the reservoir to the tubing, m^3/h
w_oil_bottom(i) = q_bottom_oil(i)*rho_oil/3600; % kg/sec

q_bottom_gas(i) = (p_bottom_pipeline_bara_gas(i) > p_bottom_riser_bara(i))*K_r_gas*sqrt((p_bottom_pipeline_bara_gas(i) - p_bottom_riser_bara(i))/G_pipeline_bottom_gas(i)); %Flow of oil from the reservoir to the tubing, m^3/h
w_gas_bottom(i) = q_bottom_gas(i)*rho_oil/3600; % kg/sec

% if p_bottom_pipeline_bara(i) > p_gas_flow_noflow(i),    
%     w_oil_bottom(i) = w_bottom(i); % kg/s - flow of oil from pipeline to riser
%     w_gas_bottom(i) = 0; % kg/s - flow of gas from pipeline to riser
% else
%     w_oil_bottom(i) = m_oil_pipeline/(m_oil_pipeline + m_gas_pipeline)*w_bottom(i); % kg/s - flow of oil from pipeline to riser
%     w_gas_bottom(i) = 0.5*(w_bottom(i) - w_oil_bottom(i)); % kg/s - flow of gas from pipeline to riser
% end

% % Flow across bottom point
% G_pipeline_bottom(i) = rho_m_pipeline(i)/rho_water;  
% 
% q_bottom(i) = (p_bottom_pipeline_bara(i) > p_bottom_riser_bara(i))*K_r*sqrt((p_bottom_pipeline_bara(i) - p_bottom_riser_bara(i))/G_pipeline_bottom(i)); %Flow of oil from the reservoir to the tubing, m^3/h
% w_bottom(i) = q_bottom(i)*rho_m_pipeline(i)/3600; % kg/sec
% 
% if p_bottom_pipeline_bara(i) > p_gas_flow_noflow(i),    
%     w_oil_bottom(i) = w_bottom(i); % kg/s - flow of oil from pipeline to riser
%     w_gas_bottom(i) = 0; % kg/s - flow of gas from pipeline to riser
% else
%     w_oil_bottom(i) = m_oil_pipeline/(m_oil_pipeline + m_gas_pipeline)*w_bottom(i); % kg/s - flow of oil from pipeline to riser
%     w_gas_bottom(i) = 0.5*(w_bottom(i) - w_oil_bottom(i)); % kg/s - flow of gas from pipeline to riser
% end

% Flow across production choke
C_pc = 50*50^(u_PC_simple(i)-1);
K_pc = 12*0.865*C_pc;
G_riser_top(i) = rho_m_riser(i)/rho_water;

q_pc(i) = (p_top_riser_bara(i) > p_sep_bara)*K_pc*sqrt((p_top_riser_bara(i) - p_sep_bara)/G_riser_top(i)); %Flow of oil from the reservoir to the tubing, m^3/h
w_pc(i) = q_pc(i)*rho_m_riser(i)/3600; % kg/sec

if p_bottom_riser_bara(i) < p_bottom_riser_onlyoil_bara(i),
    w_oil_pc(i) = m_oil_riser/(m_oil_riser + m_gas_riser)*w_pc(i); % kg/s - flow of oil from riser to separator
    w_gas_pc(i) = w_pc(i) - w_oil_pc(i); % kg/s - flow of gas from riser to separator
else
    w_oil_pc(i) = w_pc(i); % kg/s - flow of oil from riser to separator
    w_gas_pc(i) = 0; % kg/s - flow of gas from riser to separator
end

% The time derivative of the states
dm_gas_riser = w_gas_bottom(i) - w_gas_pc(i);
dm_oil_riser = w_oil_bottom(i) - w_oil_pc(i);
dm_gas_pipeline = w_gas_in(i) - w_gas_bottom(i); 
dm_oil_pipeline = w_oil_in(i) - w_oil_bottom(i);
  
dm = [dm_gas_riser; dm_oil_riser; dm_gas_pipeline; dm_oil_pipeline];    % The derivatives of the states in the system