%%%%Simple Static Model%%%
clc
clear all


g = 9.81;                           %Gravity (m/s2)
Wg_in=0.0024;                       %Inlet mass flow rate of gas (Kg/sec)
Wl_in=0.39298;                      %Inlet mass flow rate of liquid (Kg/sec)
W = Wg_in+Wl_in;                    %Inlet mass flow rate (Kg/sec)
R = 8314;                           %Gas constant (J/(K.Kmol))
M_g = 29;                           %Molecular weight of Gas (kg/kmol)

p_s=101325;                         %Seperator pressure (pa)
p_vmin = 0;                         %minimum Pressure drop over the valve (Pa)
T=15+273.15;                        %Riser temprature (K)
par.r2 = 0.025;                     %Radius of riser (m)
par.A2 = pi*par.r2^2;               %Cross section area of riser (m2)
rho_g= (p_s+p_vmin)*M_g/(R*T);      %Average gas density at the outlet(Kg/m3)
rho_l=1000;                         %Liquid density (Kg/m3)

alpha_g = Wg_in/(Wg_in+Wl_in);                                      %Average gas mass fraction
alpha_l=(1-alpha_g)*rho_g/((1-alpha_g)*rho_g+alpha_g*rho_l);        %liquid volume fraction

rho = alpha_l*rho_l+(1-alpha_l)*rho_g;        %Average mixture density
L_r = 5.15;                                   %Length of riser
z_star= 0.26;                                 %Bifurcation point
cd = 0.31;                                    %Discharge coefficient of valve
K_pc = sqrt(2)* par.A2;                       %Valve constant (m2)
a = (1/rho)*((W/K_pc)^2);                     %Constant parameter
p_star= (rho_l*g*L_r)+p_s+ p_vmin;            
fz_star = z_star*cd/sqrt(1-z_star^2*cd^2);                             
delta_p_star = a/fz_star^2;
p_fo = p_star-delta_p_star;                   %inlet pressure at fully open position of the valve(at z=1) (pa)


z_t = 0.2:0.001:1;               %Different valve openings
n = length(z_t);
Pin = zeros(1,n);                %Inlet Pressure
K_z_t = zeros(1,n);              %Static gain of the system (pa)



for i = 1:n
    fz = z_t(i)*cd/sqrt(1-z_t(i)^2*cd^2);
    Pin(i) =(a/fz^2 + p_fo)/1000;
    K_z_t(i) = (-2*a/z_t(i)^3*cd^2)/1000; 
end

figure(1)
clf
subplot(2,1,1);
plot(z_t,Pin,'k','LineWidth',2);
xlabel('Z');
ylabel('Inlet Pressure [kPa]');
hold on
grid on

subplot(2,1,2);
plot(z_t,K_z_t,'k','LineWidth',2);
xlabel('Z');
ylabel('K(z)');
hold on
grid on
