% June 2017. INFO: uses event function to stop integrating when dew point 
%is reached. 
% 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; Vliq-dummy liquid volume used to start 2P 
% 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]

% OBS. tol in event function needs to be smaller then tol used in HEX.m
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 = 15; t1 = 0; t2 = 0;
Z1P = 0.9421;                          % at 67  degC
ZG0 = 0.8519;  ZL0 = 0.06;             % at -11.1 
v0 = 1;                                % at -11.5 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.1 degC, from dew point calculation
Tdew = -11.1+273.15; Tbub = -107.9+273.15;   % at 1.8234 MPa

% load initial conditions for 1 phase (vapor) simulation 
load w01P-10Q0 % 10 cells

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 = w1(end,2*M+NC*M+1:3*M+NC*M);  % temperature (same in both phases), [K]

y0=zeros(1,NC*M); 
for j=1:M
   y0(j*NC-(NC-1):j*NC) = z0(1:4);
end

Zg0 = ones(1,M)*Z1P;               % gas compressibility
P00 = linspace(Pin,Pin-DP,M+2);    % Pressure
P0 = P00(2:end-1);x0 = zeros(1,NC*M);

for j=1:M
    x0(j*NC-(NC-1):j*NC) = z0(1:4);
end

VL0 = w1(end,3*M+3*M*NC+1:4*M+3*M*NC); 
Zl0 = ones(1,M)*0.06;
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,'Events',@efun);

%ODE solver
tstart = 0;
tfin = 200;
to=tstart;
wo=[];
teo=[];
weo=[];
ieo=[];

for i=1:2 % i = number of simulations
   % parameters passed to HEX.m 
   param = [M,NC,V,F0,Pin,DP,Tin,Qmax,t1,t2,Z1P,ZG0,ZL0,v0,z0,xin, yin, Tdew, Tbub]; 
   % integrate 
   [t,w,te,we,ie] = ode15s(@(t,w) HEX_T(t,w,param),[tstart tfin], w0, options);
   lt = length(t);

   %save output
   to = [to; t(2:lt)];  % integration time
   wo = [wo; w(2:lt,:)];% output
   teo = [teo; te];     % time when event is triggered
   weo = [weo; we];     % values of the solutions when the event is triggered
   ieo = [ieo, ie];     % number of cell where the event is triggered
         
   % conditions for 2nd simulation
   tstart = te;         % start a new simulation after the event is triggered
   LC = w(end,M:M:end);                   % results for cell # M
   A=ones(M,1);B=eye(4);C=repmat(B,1,M)';
   D=blkdiag(A,C,A,A,C,C,A,A,A,A);        % selection matrix to pass conditions from last cell of 1st simulation as initial conditions for all cells in 2nd simualtion
   w0 = D*LC';                            % new set of initial conditions
   Pin = LC(end);                         % new pressure, to keep the linear pressure drop across HEX
   %Tin=Tdew;
   Tin=LC(7);
   w0(6*M+3*M*NC+1:7*M+3*M*NC) = linspace(Pin,Pin-DP,M);
end

% extracting states from both simualtions
 N = wo(:,1:M)';                              % total holdup in [kmol]  
 n = wo(:,M+1:M+M*NC)';                       % component holdup,[kmol]
 U = wo(:,M+NC*M+1:2*M+NC*M)';                % total internal energy, [MJ]
 T = wo(:,2*M+NC*M+1:3*M+M*NC)';              % temperature        [K]
 x = wo(:,3*M+M*NC+1:3*M+2*M*NC)';            % liquid composition
 y = wo(:,3*M+2*M*NC+1:3*M+3*M*NC)';          % vapor composition
 VL = wo(:,3*M+3*M*NC+1:4*M+3*M*NC)';         % liquid volume [m3]
 Zg = wo(:,4*M+3*M*NC+1:5*M+3*M*NC)';         % compress  gas [-]
 Zl = wo(:,5*M+3*M*NC+1:6*M+3*M*NC)';         % compress liquid [-]
 P = wo(:,6*M+3*M*NC+1:7*M+3*M*NC)';          % pressure, [MPa] 

figure
plot(to(2:end),T)
grid on
xlabel('Time [min]')
ylabel('Tempertaure [K]')



