
% Fil nycobatch.m  
% Programmet simulerer en vanlig batch kolonne med nivaaregulering av 
% kondensatortanken, og sammensetningsregulering ved hjelp av temperatur.
% Programmet bruker funksjonene sumy.m og ycomp.m av Katrine Hilmen og 
% wilsonv.m av Steinar Hauan.
% Kall i matlab:
% Kjoer foerst lukket batch for aa oppnaa en profil i kolonnen
% [t,x]=ode15s('nycobatch',[0 50],A)
% Programmet er laget av Per Furu og Ketil Eik

function xdot = nycobatch(t,A)

global Xglob Tkglob

% Antall trinn inkludert koker og kondensator

NT=32;

% Antall komponenter

NC=3;



% Oevrige data
tk=Tkglob;        % Start temperatur (K)
V=121;            % Dampstroem (kmol/h)
MD=A((NC*NT)+1);  % Holdup kondensator (kmol)
MB=A((NC*NT)+2);  % Holdup koker (kmol)
MK=0.037;         % Holdup paa hvert trinn i kolonne (kmol)
MDmin=0.001;      % Brukes for aa unngaa at kondensatoren gaar tom



 
% Initsiering av sammensetning
X=A(1:NC*NT);
for i=1:NT
  x(:,:,i)=X((i-1)*NC+1:i*NC);
end


% Damp-vaeske likevekt (bruker funksjonen wilsonv.m, ycomp.m og sumy.m)
y=zeros(NC,1,NT);

Tkut=[];
for i=1:NT-1;    
Xglob=x(:,:,i);
Tk=fzero('sumy',tk(i));
Tkut(i)=Tk;
y(:,:,i)=ycomp(Tk);
Tkut=[Tkut];
end
Tkglob=Tkut;

% Regulering av sammensetningen, temperatur regulator

KC = -50;                              % Forsterkning
Tks = 337.75;                          % Settpunkt temperatur  (K)
err= Tkut(31)-Tks;                     % Proporsjonalfeil
L0I=0;                                 % Vaeskestroem (kmol/h) 
taui=0.5;                              % Integraltid
interr = A((NC*NT)+3);                 % Integralfeil
LT0 = L0I - KC*err;% - KC/taui*interr; % Regulering av vaeskestroemen
LT=LT0*(1-exp(-(MD/MDmin))) % Brukes for at kondensator ikke skal gaa tom

if LT < 0
   LT  = 0.01;
   err = 0;     % Skrur av integrasjonen for aa motvirke windup
end

% Dampstroem 
i=1:NT-1;    
V(i)=V*ones(1,NT-1);

% Vaeskestroem
i=2:NT;
L(i)=LT*ones(1,NT-1);

% Total kondensator
y(:,:,NT)=x(:,:,NT);

% Holdup i kolonnen
i=2:NT-1;
M(i)=MK*ones(1,NT-2);

% Holdup i kokeren
MBdot=L(2)-V(1);

% Nivaaregulering av holdup i kondensatortank P-regulering

KCh = 10;
MDs= 2.5; % kmol
D0= V(NT-1)-L(NT)
errh = MD-MDs;
D= D0 + KCh*errh

if D < 0
  D = 0;
  errh = 0;

end



% Holdup i kondensatoren
MDdot=V(NT-1)-L(NT)-D;

% Dynamiske ligninger (komponentbalanser)
Mxdot=zeros(NC,1,NT);
xdot =zeros(NC,1,NT);

% Kokeren
Mxdot(:,:,1)= L(2)*x(:,:,2)-V(1)*y(:,:,1);
xdot(:,:,1)=(Mxdot(:,:,1)-x(:,:,1)*MBdot)/MB;

% Kolonnen
for i=2:NT-1
xdot(:,:,i)=(L(i+1)*x(:,:,i+1)-L(i)*x(:,:,i)+V(i-1)*y(:,:,i-1)-V(i)*y(:,:,i))/M(i);
end


% Kondensatoren
Mxdot(:,:,NT) = V(NT-1)*y(:,:,NT-1)-L(NT)*x(:,:,NT)-D*x(:,:,NT);
xdot(:,:,NT)=(Mxdot(:,:,NT)-x(:,:,NT)*MDdot)/MD;

% Utverdier av holdup i koker og kondensator

MDut=MDdot;
MBut=MBdot;

%Integrator i regulatoren
errut=err;

% Utverdier av sammensetningene paa hvert trinn
xdotut=xdot(:);

% Funksjonen returnere en kolonnevektor 
(trinn 1 til NT med NC komponenter,MD,MB ,err og D)
xdot=[xdotut;MDut;MBut;errut;D]; % Ved bruk av PI-regulering

%xdot=[xdotut;MDut;MBut;D]; % ved bruk av P-regulering








