clear
% Solution to the 2001 Final Project

% This program solves the convective diffusion equation for
% the mixing of two streams in a microchannel, where the species
% undergo convection in the x-direction, diffusion in the y
% direction, and the reversible reaction A + B -- >  AB.  The
% problem and equations are laid out in the problem handout.  The
% solution approach is to discretize the domain in the y-direction
% and solve the parabolic PDE as a system of 1st order ODE's in the
% x (flow) direction using the Euler method.  The x derivatives are
% calculated in the attached function cderiv.m.

% The objective of the program is to determine the optimum initial
% width of the two inlet streams, where the diffusivity of the 
% two reactants is quite different.  Basically, we determine whether
% the high mobility stream should be narrow at a high concentration
% and the low mobility stream wide at a lower concentration, or whether
% the reverse strategy is more appropriate.

% We pass some parameters as global variables to the function which
% calculates the derivatives:
global db dab kstar hatta

% We set the discretization in the y-direction:
n=40;
dy=1/n;
y=[0:dy:1];  % These are the values of y where we compute concentrations.

% We set the dimensionless parameters:
db=.1;     % The dimensionless diffusivity of b (e.g., Db/Da)
dab=.1;    % The dimensionless diffusivity of ab
kstar=100; % The equilibrium coefficient
hatta=10;  % The hatta number (reaction/diffusion)
cb0=1;     % The dimensionless total amount of B introduced

% We are interested in the rate of approach to equilibrium.  The
% equilibrium concentration of AB is given by:

cabeq=(1+cb0)/2+1/2/kstar-(((1+cb0)/2+1/2/kstar)^2-cb0)^.5;


% We will solve for an array of values of w, the width of stream A:
wall=[.1:.1:.9];

% We are interested in doing contour plots for specified values of
% w only:
wpick=[.2 .5 .8];

figure(1)

% We will solve for each value of w in sequence:
for kk=1:length(wall)
w=wall(kk);

% OK, now we initialize the arrays and integrate:
cab=zeros(n+1,1);
ca=zeros(n+1,1);
cb=zeros(n+1,1);

% We set the initial concentrations:
i=find(y < w);  % These are the indices where y is less than w
ca(i)=ones(size(ca(i)))/w;
i=find(y > w);  % The indices where y is greater than w
cb(i)=cb0*ones(size(cb(i)))/(1-w);
i=find(y==w); % The index (if any) where y=w
% Note that if such an index occurs, both ca and cb are multivalued.
% Thus, we simply divide the concentration by two - it's a bit more
% accurate.
if i
 ca(i)=1/w/2;
 cb(i)=cb0/(1-w)/2;
end

% We also patch things up a bit by rescaling ca and cb to make sure there's
% the right amount of each in the initial stream.  This is a "fix" which is
% zero if the discretization in the y direction is high enough.  Its primarily
% important for w close to zero or one, where one of the streams is narrow.

% We simply renormalize the streams by the computed averages:
wt=ones(1,n+1);
wt(1)=0.5;
wt(n+1)=0.5;
wt=wt/n;
caave=wt*ca;
cbave=wt*cb;

ca=ca*1/caave;
cb=cb*cb0/cbave;

% OK, now we're ready to integrate forward in x.  We will use an
% explicit integration method, so we require dx to be small for
% numerical stability:

dx=0.4*dy^2;

% we want to save the values of x, ca, cb, and cab for later plotting:

x=[0:dx:.5];
m=length(x); % This is the number of nodes in the x direction

caall=zeros(m,n+1);
cball=zeros(m,n+1);
caball=zeros(m,n+1);

% We initialize the first row of caall and cball:
caall(1,:)=ca';
cball(1,:)=cb';

% We are interested in the approach to equilibrium, so we will track
% the ratio of the average concentration to cabeq.  This is also
% initialized:
cabratio=zeros(m,1);

% And now we perform the integration:
for i=2:m
 [cadot cbdot cabdot]=cderiv(ca,cb,cab);
 ca=ca+cadot*dx;
 cb=cb+cbdot*dx;
 cab=cab+cabdot*dx;
 
 % And we calculate the approach to equilibrium:
 cabratio(i)=wt*cab/cabeq;
 
 % And we save the concentration arrays:
 caall(i,:)=ca';
 cball(i,:)=cb';
 caball(i,:)=cab';

 % The following bit should be commented out if you don't want to see
 % the movie of the concentrations, as it slows things way down!  Set
 % flag = 1 if you want to see the movie, flag=[ ] if you don't.
 flag=[ ];
 if flag
  figure(1)
  plot(y,ca,y,cb,y,cab,y,cabratio(i)*ones(size(y)),'--')
  text(.2,.5,['x = ',num2str(x(i))])
  title('Concentration evolution in x as a function of y')
  xlabel('y')
  ylabel('dimensionless concentration')
  legend('Ca','Cb','Cab',' < Cab > /Cabeq',-1)
  drawnow
 end
end

% We want to save some bits of information for later plotting about
% the behavior for a particular value of w:

cabratioall(:,kk)=cabratio;

% And we find the value of x for which we approach 80%  of equilibrium:
[junk,i]=min(abs(cabratio-0.8));
xeq(kk)=x(i);

% For specific values of w we wish to do contour plots:
j=find(wpick==w); % This returns 1, 2, or 3 if w is one of the picked values

if j

 % We shall use nc contour lines:
 nc=10;
 
 % Now we do the contour plots:
 figure(3*j-1) % The first of three figures for this w
 [c,h]=contour(x,y,caall',nc);
 clabel(c,h)
 title(['Concentration Contours for Ca at w = ',num2str(w)])
 xlabel('x')
 ylabel('y')
 
 figure(3*j) % The second of three figures for this w
 [c,h]=contour(x,y,cball',nc);
 clabel(c,h)
 title(['Concentration Contours for Cb at w = ',num2str(w)])
 xlabel('x')
 ylabel('y')

 figure(3*j+1) % The third of three figures for this w
 [c,h]=contour(x,y,caball',nc);
 clabel(c,h)
 title(['Concentration Contours for Cab at w = ',num2str(w)])
 xlabel('x')
 ylabel('y')
end

end

% We are interested in the ratios at the three specified values of
% w given by wpick:
linetypes=str2mat('-',':','--','-.','*'); % We differentiate among the linetypes
figure
p=find(wpick(1)==wall);
plot(x,cabratioall(:,p),linetypes(1,:))
labels=['w = ',num2str(wall(p))];
hold on
 for i=2:length(wpick)
 p=find(wpick(i)==wall);
 plot(x,cabratioall(:,p),linetypes(i,:))
 labels=str2mat(labels,['w = ',num2str(wall(p))]);
end
legend(labels)
xlabel('x')
ylabel(' < cab >  /  < cab > eq')
title('Comparison of approach to equilibrium for several values of w')
hold off

% And we have the final figure:
figure
plot(wall,xeq)
title('Axial distance to approach 80%  of equlibrium')
ylabel('Distance')
xlabel('Width of Stream A')
axis([0 1 .1 .4])

disp('       w       xeq');disp([wall',xeq']);

echo on

% Based on these calculations, we find that the approach to equilibrium is
% fastest if the width of the more mobile stream is as small as possible.
% The improvement is relatively modest, with an equlibration distance of 0.2
% for w=0.1 and 0.272 at w=0.5, a 26%  decrease.  The maximum value of the 
% equilibration distance occurred at w=0.8, with a value of 0.296.  Examination 
% of the contour plots indicates another advantage for small w: the concentration 
% profile of AB is much more uniform.  This would be important if there were forms 
% of axial dispersion (unlike what was assumed here), or if there are other non-
% linear reaction effects.  For both these reasons, therefore, it is desirable
% to make the mobile stream narrow and concentrated, and the lower mobility
% stream wider.
echo off
 
 


_______________________________________


function [cadot,cbdot,cabdot]=cderiv(ca,cb,cab)
% This function takes in the three concentrations ca, cb, and cab
% and returns the time derivatives of the concentrations.  The
% concentrations are arrays, as they are discretized in space.
% They are assumed to be column vectors.

% We also bring in a number of dimensionless groups as global 
% parameters:
global db dab kstar hatta

% db is the dimensionless diffusivity ratio Db/Da
% dab is the dimensionless diffusivity ratio Dab/Da
% kstar is the equilibrium constant kf < Ca0 > /kr
% hatta is the hatta number ha=[kf < Ca0 > b^2/Da]^.5

% we calculate n, the spatial discretization, from the length of ca:
n=length(ca)-1;

% OK, now for the derivatives:

rstar=-ca.*cb+cab/kstar;

% We produce the spatial discretization matrix:
a=-2*diag(ones(n+1,1));
a=a+diag(ones(n,1),1);
a=a+diag(ones(n,1),-1);

a=a*n^2;  % Note that dy = 1/n;

cadot=a*ca+hatta^2*rstar;
cbdot=db*a*cb+hatta^2*rstar;
cabdot=dab*a*cab-hatta^2*rstar;

% These derivatives are only valid for the interior points, e.g.,
% i = 2:n.

% Now we have to fix the first and last elements of cadot, cbdot,
% and cabdot to account for the boundary conditions.  We could do
% this after the integration step, or (if the initial conditions
% satisfy the B.C.'s) put it here in the derivative function.  This
% is a superior method if you are using a higher order integration
% method in x, but doesn't matter if you are using the Euler method.
% We'll put it here to make the main program a little shorter.

% The boundary conditions are derived from zero normal derivative
% conditions.  We simply set an order dy^2 finite difference
% approximation for the derivative equal to zero and solve for ca(1)
% and ca(n+1), etc., or the x derivative of this quantity:

cadot(1)=4/3*cadot(2)-1/3*cadot(3);
cadot(n+1)=4/3*cadot(n)-1/3*cadot(n-1);

cbdot(1)=4/3*cbdot(2)-1/3*cbdot(3);
cbdot(n+1)=4/3*cbdot(n)-1/3*cbdot(n-1);

cabdot(1)=4/3*cabdot(2)-1/3*cabdot(3);
cabdot(n+1)=4/3*cabdot(n)-1/3*cabdot(n-1);

% and that's it!