% June 2017. Main code for runing a simulation for vapor phase.

% M-nr cells; NC-number of components (1)Nitrogen - (2)Methane - (3)Ethane
% - (4)Propane; V-HEX Volume,[m3]; F0-inlet flow,[kmol/min]; 
% Pin-inlet pressure,[MPa];% DP-pressure drop along HEX,[MPa]; 
% Tin-inlet temp,[K]; Qmax-heat removed,[MJ/min]; t1-time
% when Q!=0,[min]; t2-time when Q reaches max,[min]; Z1P-Z for 1P; ZG0-Z
% for gas for 2P; z0-overall molar fraction; xin-liq inlet comp for 2 phase;
% yin-vap inlet comp for 2 phase; Tdew-dew point; Tbub-buble point;
% param = [N,NC,V,F0,Pin,DP,Tin,Qmax,t1,t2,Z1P,ZG0,ZL0,v0,z0,xin, yin, Tdew,Tbub]

clear;close all;clc

M  = 10; NC = 4; V = 100; F0 = 0.0445*60; 
Pin = 18*0.101325; DP = 5*1e-3*0.1; Tin = 67+273.15;
Qmax = 10; t1 = 0; t2 = 0;
Z1P = 0.9421;                          % at 67  degC
ZG0 = 0.8527; ZL0 = 0.0674; v0 = 1;    % at -11.1 degC
z0 = [0.06 0.4 0.4 0.14]; 
xin = [0.0035    0.0656    0.3985    0.5324]; 
yin =[0.06 0.4 0.4 0.14];   % at -11.5 degC, from dew point calculation
Tdew = -11.1+273.15; Tbub = -107.9+273.15;   % mean from Pin to Pin-DP

param = [M,NC,V,F0,Pin,DP,Tin,Qmax,t1,t2,Z1P,ZG0,ZL0,v0,z0,xin, yin, Tdew, Tbub];

% load initial conditions for 1 phase (vapor) simulation 

load w01P-10Q0 % 10 cells
%load w01P-60Q0 % 60 cells

Pout = Pin-DP;

nT0 = w1(end,1:M);       % total holdup in  (both phases),[kmol]  
n0 = zeros(1,NC*M);
for j=1:M                % component 1 holdup in tank (both phases),[kmol]
    n0(NC*j-3:NC*j) = nT0(j).*z0;
end

U0 = w1(end,M+NC*M+1:2*M+NC*M);  % total internal energy (both phases),[MJ]
T0 = ones(1,M)*Tin;              % temperature (same in both phases), [K]

y0=zeros(1,NC*M);%vapor composition = overall compositions(no phase change)
for j=1:M
   y0(j*NC-(NC-1):j*NC) = yin(1:4);
end

Zg0 = ones(1,M)*Z1P;            % gas compressibility
P00 = linspace(Pin,Pin-DP,M+2); % pressure inlet to outlet, [MPa]
P0 = P00(2:end-1);              % pressure inside the HEX

x0 = zeros(1,NC*M);             % dummy liquid compositions
for j=1:M
    x0(j*NC-(NC-1):j*NC) = xin(1:4);
end

VL0 = w1(end,3*M+3*M*NC+1:4*M+3*M*NC); % dummy liquid volume
Zl0 = ones(1,M)*0.06;                  % dummy liquid comprsibility  

w0 = [nT0,n0,U0,T0,x0,y0,VL0,Zg0,Zl0,P0]; % vector of initial condition

%ODE solver options********************************************************

nrDyn = NC+2;          %number of dynamic equations per cell
nrAlg = 2*NC+5;        %number of algebraic eq per cell
d1 = ones(1,nrDyn*M);
d2 = zeros(1,nrAlg*M);
Mass = diag([d1,d2]);
options = odeset('mass',Mass,'Stats','on');

%ODE solver;integrate with algortihm based on comparison with the dew point
[t,w]=ode15s(@(t,w) HEX_T(t,w,param),[0 200],w0,options);       


