% programa 5.13, W.L.Luyben, pag. 154
% "batch reactor simulation"
% modif. E.Quintero.M 26may'98
%-------------------------------------------------
% datos: delta=paso de integración
%        tf=tiempo final
%        tprint=intervalo de impresión
%----------------------------------------------
clear
tf=149; delta=0.002; time=0; tprint=1; 
alpha1=729.5488; alpha2=6567.587; 
avp=-8744.4; bvp=15.70036;
kc=2; ca=0.8; cb=0;
t=80; tm=80; tj=259; pj=34.4; pset=12.6;
start=1; fwo=0; full=-1; vj=0; vjtj=0;
tflag=tprint; cam=-2; ramp=0.005; 
dens=18*pj*144/(1545*(tj+460));
xs=0; xw=80;
while(time<=tf)
   while(time<=tflag)
      k1=alpha1*exp(-15000/(1.99*(t+460)));
      k2=alpha2*exp(-20000/((t+460)*1.99));
      ptt=3+(t-50)*12/200;   % transmitter
      pc=7+kc*(pset-ptt);    % controller
      if(pc>15) pc=15; end
      if(pc< 3) pc= 3; end
      xs=(pc-9)/6;           % valves
      xw=(9-pc)/6;
      if(xs>1) xs=1; end
      if(xs<0) xs=0; end
      if(xw>1) xw=1; end
      if(xw<0) xw=0; end
      if(start>=0)      % test for steam
         if(pj<35)
            ws=xs*112*sqrt(35-pj);
         else
            ws=0;
         end
         qj=-1000*56.5*(tj-tm)/60;
         wc=-qj/939;
         dendot=(ws-wc)/18.83;
         dens=dens+delta*dendot;        
         flagp=-1; flagm=-1; dtj=1;
         loop=0;
         while(loop<50)   %  iterative loop
            pj=exp(bvp+avp/(tj+460));
            loop=loop+1;
            if(loop>=50) error('steam temp loop'); end
            dcalc=18*pj*144/(1545*(tj+460));
            if(abs(dens-dcalc)<0.0011) break; end
            if(dens<=dcalc)
               if(flagm>=0) dtj=dtj/2; end
               tj=tj-dtj;
               flagp=1;
            else
               if(flagp>=0) dtj=dtj/2; end
               tj=tj+dtj;
               flagm=1;
            end
         end
      else
         fwo=100*sqrt(20)*8.33*xw/62.3;
         ws=0; wc=0; pj=0; xs=0;
         if(full<=0)
            ao=vj*56.5/18.83;
            vj=vj+delta*fwo;
            if(vj>=18.83) full=1; end
            qj=400*ao*(tm-tj)/60;
            vjtj=vjtj+delta*(fwo*80+qj);
            if(vj>0)
               tj=vjtj/vj;
            else
               tj=80;
            end
         else      % full jacket
            qj=400*56.5*(tm-tj)/60;
            tprint=2;
            delta=0.05; % use bigger delta once jacket is full of water
            tjdot=fwo*(80-tj)/18.83+qj/(18.83*62.3);
            tj=tj+delta*tjdot;
         end
      end
      cadot=-k1*ca;
      cbdot=k1*ca-k2*cb;
      qm=160*56.5*(t-tm)/60;
      tdot=(k1*ca*40000+k2*cb*50000)/50 - qm/(42.4*50);
      tmdot=(qm-qj)/(512*0.12*9.42);
      if(time==0) xx=[time ca cb xs xw t tm tj fwo qj qm pc]; end
      time=time+delta;    % integration
      ca=ca+cadot*delta;
      cb=cb+cbdot*delta;
      t=t+tdot*delta;
      tm=tm+tmdot*delta;
      if(t>300) error('temp > 300 °R'); end
      if(t>200) start=-1; cam=1; end
      if(cam>0) pset=pset-delta*ramp; end
   end
   tflag=tflag+tprint;
   xx=[xx; time ca cb xs xw t tm tj fwo qj qm pc];
end
format short e
subplot(4,2,1),plot(xx(:,1),xx(:,2)),ylabel('Ca'),grid
subplot(4,2,3),plot(xx(:,1),xx(:,3)),ylabel('Cb'),grid
subplot(4,2,5),plot(xx(:,1),xx(:,6)),ylabel('T'),grid
subplot(4,2,7),plot(xx(:,1),xx(:,7)),ylabel('Tm'),xlabel('Time, min'),grid
subplot(4,2,2),plot(xx(:,1),xx(:,8)),ylabel('Tj'),grid
subplot(4,2,4),plot(xx(:,1),xx(:,9)),ylabel('Fwo'),grid
subplot(4,2,6),plot(xx(:,1),xx(:,12)),ylabel('Pc'),grid
subplot(4,2,8),plot(xx(:,1),xx(:,11)),ylabel('Qm'),xlabel('Time, min'),grid
