% Dynamic analysis of a plug flow reactor
%
% The following equations model a plug flow reactor (PFR)
%
%   Material balance
%
%   ca  = -v*ca  + D*(ca   + (1/r)ca ) - k*ca^2                 (1)
%     t        z        rr          r
%
%   Energy balance
%
%   T  = -v*T  + k*/(rho*Cp)*(T   + (1/r)*T )                   (2)
%    t       z                 rr          r
%
%        - (dH*k/(rho*Cp)*ca^2
%
%   kr = k0*exp(-E/(R*T))                                       (3)
%
% The variables and parameters for this model are (in cgs units)
%
%   Reactant concentration          ca        (eq. (1))
%
%   Temperature                     T         (eq. (3))
%
%   Reaction rate constant          kr        (eq. (3))
%
%   Time                            t
%
%   Radial position                 r
%
%   Axial position                  z
%
%   Entering concentration          ca0          0.01
%
%   Entering temperature            T0           305
%
%   Wall temperature                TW           305
%                                                355
%
%   Reactor radius                  r0           2
%
%   Reactor length                  zl           100
%
%   Linear velocity                 v            1
%
%   Mass diffusivity                D            0.1
%
%   Thermal diffusivity     alpha = k/(rho*Cp)   0.1
%
%   Liquid density                  rho          1.0
%
%   Liquid specific heat            Cp           0.5
%
%   Heat of reaction                dH          -10,000
%                                                    
%   Specific rate constant          k0           1.5e+09
%                                                2.0e+09
%
%   Activation energy               E            15,000
%
%   Gas constant                    R            1.987
%
% Note that there are two cases, corresponding to the two
% values of the specific rate constant, k0.  This is a
% sensitive parameter since it multiplies a stongly nonlinear
% term (Arrhenius temperature dependency) in eq. (3).  Thus,
% two values are programmed k0 = 1.5e+09, 2.0e+09, to indicate
% the effect of k0.
%
% The solution to this system, ca(r,z,t) (from eq. (1)) and
% T(r,z,t) (from eq. (2)) is computed by a fixed step Euler
% integration in t.  Also, the choice of an Euler step is
% sensitive because of the nonlinearity of eq. (3).
%
% Open output file
  fid1=fopen('reactor.out','w');
%
% Model parameters (defined in the comments above)
  ca0=0.01;
  T0=305.0;
%
% Case 1: Cooled reactor wall
% Tw=305.0;
%
% Case 2: Heated reactor wall
  Tw=355.0;
  r0=2.0;
  zl=100.0;
  v=1.0;
  D=0.1;
  alpha=0.1;
  rho=1.0;
  Cp=0.5;
  dH=-10000.0;
  E=15000.0;
  R=1.987;
%
% Grid in axial direction
  nz=21;
  dz=zl/(nz-1);
  for j=1:nz
    z(j)=(j-1)*dz;
  end
%
% Grid in radial direction
  nr=5;
  dr=r0/(nr-1);
  for i=1:nr
    r(i)=(i-1)*dr;
  end
  drs=dr^2;
%
% Case 1: Low reaction rate constant
% rk0=1.5e+09;
%
% Case 2: High reaction rate constant
  rk0=2.0e+09;
%
% Series of solutions for different Euler steps
  for ncase=1:1
%
% Initial, final times
  t=0.0; tf=300.0;
%
% Integration step, number of Euler steps for each output
  if ncase==1 h=0.1; nout=100; end
%
% Initial conditions
  for j=1:nz
  for i=1:nr
    ca(i,j)=0.0;
     T(i,j)=T0;
  end
  end
%
% Write h
  fprintf(fid1,'\n\n ncase = %5d    h = %10.3e\n\n',ncase,h);
%
% Write heading
  fprintf(fid1,'     t   ca(r,zl,t)  T(r,zl,t)\n');
%
% Integrate until t = tf
  while t < tf*1.0001
%
% Monitor solution by displaying t
  t
%
% Write selected output (a quadrature, e.g., Simpson's
% rule, could be added here for the average output ca, T)
  for i=1:nr
    fprintf(fid1,'%6.1f%12.6f%10.2f\n',t,ca(i,nz),T(i,nz));
  end
  fprintf(fid1,'\n');
%
% Take nout Euler steps
  for iout=1:nout
%
%   Temporal derivatives
%
%     Cover the nr x nz grid
%
%       Entering conditions (z = 0)
        for i=1:nr
          ca(i,1)=ca0;
          cat(i,1)=0.0;
          T(i,1)=T0;
          Tt(i,1)=0.0;
        end
%
%       Rest of reactor
        for j=2:nz
%
%         Centerline (r = 0)
          rk=rk0*exp(-E/(R*T(1,j)));
          cat(1,j)=-v*(ca(1,j)-ca(1,j-1))/dz...
                  +4.0*D*(ca(2,j)-ca(1,j))/drs...
                  -rk*ca(1,j)^2;
          Tt(1,j)=-v*(T(1,j)-T(1,j-1))/dz...
                  +4.0*alpha*(T(2,j)-T(1,j))/drs...
                  -dH*rk/(rho*Cp)*ca(1,j)^2;
%
%         Wall (r = r0)
          rk=rk0*exp(-E/(R*Tw));
          cat(nr,j)=-v*(ca(nr,j)-ca(nr,j-1))/dz...
                   +2.0*D*(ca(nr-1,j)-ca(nr,j))/drs...
                   -rk*ca(nr,j)^2;
           T(nr,j)=Tw;
          Tt(nr,j)=0.0;
%
%         Interior, r ~= 0 and r ~= r0
          for i=2:nr-1
            rk=rk0*exp(-E/(R*T(i,j)));
            cat(i,j)=-v*(ca(i,j)-ca(i,j-1))/dz...
                   +D*(ca(i+1,j)-2.0*ca(i,j)+ca(i-1,j))/drs...
                   +D*(1.0/r(i)*(ca(i+1)-ca(i-1)))/(2.0*dr)...
                   -rk*ca(i,j)^2;
            Tt(i,j)=-v*(T(i,j)-T(i,j-1))/dz...
                   +alpha*(T(i+1,j)-2.0*T(i,j)+T(i-1,j))/drs...
                   +alpha*(1.0/r(i)*(T(i+1,j)-T(i-1,j)))/(2.0*dr)...
                   -dH*rk/(rho*Cp)*ca(i,j)^2;
%
%         Next r
          end
%
%       Next z
        end
%
%   All temporal derivatives are computed
%
%   Take Euler step
    for j=1:nz
    for i=1:nr
      ca(i,j)=ca(i,j)+cat(i,j)*h;
       T(i,j)= T(i,j)+ Tt(i,j)*h;
    end
    end
    t=t+h;
%
% Next Euler step
  end
%
% Next output
  end
%
% Next case
  end