
% Fil nyco7.m  
% Programmet simulerer en lukket batch kolonne med variabel holdup i koker og
% kondensatortank inderekte regulert av sammensetningen som blir regulert
% med en temperatur regulator (P- eller PI-regulator)
% Programmet bruker funksjonen sumy.m, ycomp.m av Katrine Hilmen og wilsonv.m
% av Steinar Hauan 
% Kall i matlab:
% Kjoer foerst en av initfilen, se de forskjellige initfilen for de ulike
% stroemmene
% [t,x]=ode15s('nyco7',[0 50],A)
% Programmet er laget av Per Furu og Ketil Eik

function xdot = nyco7(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 kondensatortank (kmol)
MB=A((NC*NT)+2);  % Holdup koker (kmol)
MK=0.037;         % Holdup paa hvert trinn i kolonne (kmol)
MDmin=1;          % Minimumskarv til mengde i kondensatortank


 
% 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 = -70;                             % Forsterkning
Tks = 337.7;                          % Settpunkt temperatur (K)
err= Tkut(16)-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.001;
   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);

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

% 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);
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 og err)
xdot=[xdotut;MDut;MBut;errut]; % ved bruk av PI-regulator

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








