% Jacketed exothermic CSTR (Luyben's CSTR)
% This program simulates the noisothermal CSTR
% Disturbance is step change in feed composition
%
%initialization
  ttime=0;trand_dat=[];
% Design parameters
 Volume=48.;
%initial state condition
 density=50;cp=.75;
 Ca=0.245;T=600;Tset=T;Vset=48.;Ca0=.5;F0=40;% Reactor
 Tj0=530;Tj=594.64;T0=530;Fj0=49.9;Ucoef=150;Area=250;% Jacket
 V_jacket=3.85;dens_j=62.3;
 ko=7.08e+10;Energy=30000;Rcoef=1.99;lamda=-30000;%Reaction
 Kc1=5.0; tauT=.8; Kc2=10.;tauV=10;errT=0;errV=0;errintT=0;errintV=0;%controller gains
% Parameter setting for Euler method;
 delta_t=.01; endtime=10;samp_t=.1;
% Disturbance occurs
 Ca0=.55
%
  VCa=Volume*Ca;VT=Volume*T;
  trand_dat=[ttime Ca T Volume F0 Tj Fj0]
  tsamp_t=samp_t;
  %Simulation Starts 
  while ttime<endtime,
  % Controller
   errT=Tset-T;
   Fj=Fj0-Kc1*(errT+errintT/tauT);
  % Fj=Fj0;
   errV=Vset-Volume;
  F=F0-Kc2*(errV+errintV/tauV);
  % F=F0;
  %
  Reac_rate=ko*exp(-Energy/(Rcoef*T));
  Q=Ucoef*Area*(T-Tj);
  % Evaluation of derivative
   Volume_dot=F0-F;
   VT_dot=F0*T0-F*T-(lamda*Volume*Reac_rate*Ca+Q)/(density*cp);
   VCa_dot=F0*Ca0-F*Ca-Volume*Reac_rate*Ca;
   Tj_dot=Fj*(Tj0-Tj)/V_jacket + Q/(V_jacket*dens_j);
   % Euler method
  ttime=ttime+delta_t;
  Volume=Volume+Volume_dot*delta_t;
  VT=VT+VT_dot*delta_t;
  VCa=VCa+VCa_dot*delta_t;
  Tj=Tj+Tj_dot*delta_t;
  Ca=VCa/Volume;
  T=VT/Volume;
  errintT=errintT+errT*delta_t;
  errintV=errintV+errV*delta_t;
 if ttime>=tsamp_t,
   xtemp=[ttime Ca T Volume F Tj Fj]
   trand_dat=[trand_dat;xtemp];
   tsamp_t=tsamp_t+samp_t;
   end;
  end;
 % this is the end of simulation
  subplot(211);plot(trand_dat(:,1),trand_dat(:,2));
  title('concentration in the reactor');
  subplot(212);plot(trand_dat(:,1),trand_dat(:,3));
  title('reactor temp');
  end;
    