% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt%% This program is free software; you can redistribute it and/or% modify it under the terms of the GNU General Public License as% published by the Free Software Foundation; either version 2, or (at% your option) any later version.%% This program is distributed in the hope that it will be useful, but% WITHOUT ANY WARRANTY; without even the implied warranty of% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU% General Public License for more details.%% You should have received a copy of the GNU General Public License% along with this program; see the file COPYING. If not, write to% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,% MA 02111-1307, USA.globalepsbaNafinNtfinQfDak1ABintAintRintnptsDacafxaxa=0.75;Rg=8.314;% J/K molP=4*1.013e5;% N/m^2T=550;% Kctf=P/(Rg*T)*1e-6;% mol/cm^3Rp=0.45;% cm radius of catalyst particlea=Rp/3;Nafin=10;% mol/secNIf=10;Ntf=Nafin+NIf;caf=ctf*0.5;Qf=Nafin/caf;k1=2.25e5;%cm^3/mol sDa=0.008;% cm^2/srhob=0.60;% g/cm^3 bed densityrhop=0.68;% g/cm^3 particle densityepsb=1-rhob/rhop;% bed porosity, dimensionless%% global collocation%npts=25;[RABQ]=colloc(npts-2,'left','right');R=R*Rp;A=A/Rp;B=B/(Rp*Rp);Q=Q*Rp;Aint=A(2:npts-1,:);Bint=B(2:npts-1,:);Rint=R(2:npts-1);%% find the pellet profile at bed inlet%ca0=logspace(log10(caf)-2,log10(caf),npts)';tol=1e-12;opts=optimset('TolFun',tol);[ca,fval,info]=fsolve('pellet',ca0,opts);info;%% march down the bed%nvs=100;vfinal=4e5;vsteps=linspace(0,vfinal,nvs)';vout=vsteps;y0=[Nafin;ca];ydot0=zeros(length(y0),1);res=bed(y0,ydot0,0);ydot0(1)=-res(1);ymin=min(y0);opts=odeset('AbsTol',sqrt(eps),'RelTol',1e-10,'Events',@stop);[vout,y]=ode15i(@bed,vsteps,y0,ydot0,opts);nout=length(vout);if(nout==nvs)fprintf('hey, did not reach final conversion, increase vfinal\n');endxa=(Nafin-y(end,1))/Nafin;Naf=y(:,1);vplot=vout/1000.;%litVR=vplot(end);tableex=[vplot,Naf];%% pick out some good length locations for pellet profiles%rowsc=[1,nout];colsc=[2:npts+1];ca=y(rowsc,colsc)';table2=[R,ca];Naout=(1-xa)*Nafin;Natop=(1-xa+0.10)*Nafin;% dashedlines = [0, Naout, VR, 0, 300, Naout, VR, 2 ;% 1.1*VR, Naout, VR, Natop, 400, Naout, VR, 3 ];% Compare to the two approximations given in ch7,% Example 7.5, Figure 7.26par.k=k1;par.Nafin=Nafin;par.T=T;par.rhop=rhop;par.rhob=rhob;par.Da=Da;par.Rp=Rp;par.Rg=82.06;par.P=4;par.n=2;par.xa=xa;par.nvs=nvs;par.vfinal=vfinal;% solve reactor with: eta = 1./Phi*( 1./tanh(3*Phi) - 1/(3*Phi) );par.eta=(@(x)1./x*(1./tanh(3*x)-1/(3*x)));[vap1,xap1]=pbrsolve(par);% solve reactor with: eta = 1./Phi;par.eta=(@(x)1./x);[vap2,xap2]=pbrsolve(par);vap1=vap1/1000.;VRap1=vap1(end);vap2=vap2/1000.;VRap2=vap2(end);tableap1=[vap1xap1];tableap2=[vap2xap2];dashedlines=...[0,Naout,VR,0,VRap1,0,VRap2,0;...1.1*VR,Naout,VR,Natop,VRap1,Natop,VRap2,Natop];savefb2colloc.dattableextableap1tableap2dashedlinesif(~strcmp(getenv('OMIT_PLOTS'),'true'))% PLOTTING%plot the molar flow versus reactor volume and 75% conversion lineplot(vplot,Naf,...vap1,xap1,...vap2,xap2,...dashedlines(:,1),dashedlines(:,2),...dashedlines(:,3),dashedlines(:,4),...dashedlines(:,5),dashedlines(:,6),...dashedlines(:,7),dashedlines(:,8))axis([0,400,2,10])% TITLEend% PLOTTING