% filename run_radialOC.m
% This is the rhs of the PDEs for a reactor with radial dispersion.
% They are solved with orthogonal collocation in cylindrical geometry.

global alphap betap gammap Pem Peh Damk Biot tsurr
global xxcoll aacoll bbcoll qinvcoll wwcoll ncol

% set ncol before calling
% set parameters
alphap = 1.
betap = 0.6666667
gammap = 20.
Pem = 1
Peh = 1
Damk = 0.3
Biot = 20 %1
tsurr = 1. %0.92

% set the time interval for integration
zspan = [0 1]
options=odeset('RelTol',1.e-12,'AbsTol',1.e-6)

% call the collocation matrices, for n interior points, cylindrical geometry
%ncol = 1
ncol
coll(ncol,2,0);

% set the initial conditions for all variables
for k = 1:ncol
yi(1,k) = 1.;
yi(2,k) = 1.;
end

% move variables to one long list
num = 0;
for nvar = 1:2
for k=1:ncol
num=num+1;
y0(num) = yi(nvar,k);
end
end
y0

[z,solution] = ode45('radialrhsOC',zspan, y0,options);

% move the variables into variable names
num = 0;
for nvar = 1:2
for k=1:ncol
num=num+1;
ytemp(:,nvar,k) = solution(:,num);
end
end
for k=1:ncol
conc(:,k) = ytemp(:,1,k);
temp(:,k) = ytemp(:,2,k);
end
nn = size(conc);
for kkk = 1:nn(1)
for k=1:ncol
c(k) = conc(kkk,k);
t(k) = temp(kkk,k);
end
concb = boundvalue(c,0,0);
tempb = boundvalue(t,Biot,tsurr);
sum1 = concb*wwcoll(ncol+1);
sum2 = tempb*wwcoll(ncol+1);
for k=1:ncol
sum1 = sum1 + wwcoll(k)*c(k);
sum2 = sum2 + wwcoll(k)*t(k);
end
cavg(kkk) = 2.*sum1;
tavg(kkk) = 2.*sum2;
end
plot(z,cavg,'-r',z,tavg,'-b')
xlabel('axial location')
legend('average concentration','average temperature')
title('Problem One')


*******************************************************
% filename radialrhsOC.m
% This is the rhs of the PDEs for a reactor with radial dispersion.
% They are solved with orthogonal collocation in cylindrical geometry.

function ydot=radialrhs(t,y)
global alphap betap gammap Pem Peh Damk Biot tsurr
global xxcoll aacoll bbcoll qinvcoll wwcoll ncol

% move y into variable names

num = 0;
for nvar = 1:2
for k=1:ncol
num=num+1;
ytemp(nvar,k) = y(num);
end
end
conc = ytemp(1,:);
temp = ytemp(2,:);
concb = boundvalue(conc,0.,0.);
tempb = boundvalue(temp,Biot,tsurr);

% now compute the right-hand sides, for the time derivatives

% first do the transport terms
for j = 1:ncol
sum1 = bbcoll(j,ncol+1)*concb;
sum2 = bbcoll(j,ncol+1)*tempb;
for i=1:ncol
sum1 = sum1 + bbcoll(j,i)*conc(i);
sum2 = sum2 + bbcoll(j,i)*temp(i);
end
ydot(j) = (alphap/Pem)*sum1;
ydot(j+ncol) = (alphap/Peh)*sum2;
end

% next do the reaction terms

for j=1:ncol
rate1 = Damk*conc(j)*exp(gammap*(1.-1./temp(j)));
ydot(j) = ydot(j) - rate1;
ydot(j+ncol) = ydot(j+ncol) + betap*rate1;
end

ydot = ydot';

*******************************************************
% boundvalue.m
% program to find the value at the N+1-th collocation point
% when the values at the other collocation points are known.
% BB is the Biot number for heat, zero for mass

function yb = boundvalue(y,BB,tsurr)
global xxcoll aacoll bbcoll qinvcoll wwcoll ncol

sum = 0.;
for i=1:ncol
sum = sum + aacoll(ncol+1,i)*y(i);
end
yb = (BB*tsurr - sum)/(aacoll(ncol+1,ncol+1) + BB);


*******************************************************
% plot_radialOC.m
% plots results of the reactor calculation
global xxcoll aacoll bbcoll qinvcoll wwcoll ncol
global alphap betap gammap Pem Peh Damk Biot tsurr

ncol=6
% set parameters
alphap = 1.
betap = 0.6666667
gammap = 20.
Pem = 1
Peh = 1
Damk = 0.3
Biot = 1
tsurr = 0.92

% run for six z
dataz = [0.2 0.4 0.5 0.6 0.7 1.0]
% set the time interval for integration
zspan = [0 dataz(1)]

% plotting points
x = 0:.025:1;

% call the collocation matrices, for n interior points, cylindrical geometry
coll(ncol,2,0);

for i=1:ncol
yi(1,i) = 1.;
yi(2,i) = 1.;;
end

% move variables to one long list
num = 0;
for nvar = 1:2
for k=1:ncol
num=num+1;
y0(num) = yi(nvar,k);
end
end

for kkk=1:6
[z,solution] = ode45('radialrhsOC',zspan, y0);

% move the variables into variable names
num = 0;
for nvar = 1:2
for k=1:ncol
num=num+1;
ytemp(:,nvar,k) = solution(:,num);
end
end
nn = size(ytemp)
for i=1:ncol
c(kkk,i) = ytemp(nn(1),1,i);
t(kkk,i) = ytemp(nn(1),2,i);
end
conc = c(kkk,:);
temp = t(kkk,:);
concb = boundvalue(conc,0.,0.);
tempb = boundvalue(temp,Biot,tsurr);
c(kkk,ncol+1)=concb;
t(kkk,ncol+1)=tempb;
conc(ncol+1) = concb;
temp(ncol+1) = tempb;
dd = qinvcoll*conc';
cinterp(kkk,:) = coll_interp(ncol,dd,x);
dd = qinvcoll*temp';
tinterp(kkk,:) = coll_interp(ncol,dd,x);
% set parameters for next interval
if kkk==6, break, end
zspan = [dataz(kkk) dataz(kkk+1)]
for i=1:2*ncol
y0(i) = solution(nn(1),i);
end
clear ytemp solution conc temp
end

% plot the solutions
plot(xxcoll,c(1,:),'ob',xxcoll,c(2,:),'og',xxcoll,c(3,:),'or',xxcoll,c(4,:),'oc',xxcoll,c(5,:),'om',xxcoll.c(6,:),'ok',x,cinterp(1,:),'-b',x,cinterp(2,:),'-g',x,cinterp(3,:),'-r',x,cinterp(4,:),'-c',x,cinterp(5,:),'m',x,cinterp(6,:),'-k')
xlabel('radial location')
ylabel('concentration')
legend('z = 0.2','z = 0.4','z = 0.5','z = 0.6','z = 0.7','z = 1.0')
title('Problem One')
pause
plot(xxcoll,t(1,:),'ob',xxcoll,t(2,:),'og',xxcoll,t(3,:),'or',xxcoll,t(4,:),'oc',xxcoll,t(5,:),'om',xxcoll,t(6,:),'ok',x,tinterp(1,:),'-b',x,tinterp(2,:),'-g',x,tinterp(3,:),'-r',x,tinterp(4,:),'-c',x,tinterp(5,:),'m',x,tinterp(6,:),'-k')
xlabel('radial location')
ylabel('temperature')
legend('z = 0.2','z = 0.4','z = 0.5','z = 0.6','z = 0.7','z = 1.0')
title('Problem One')