% 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.globalk_mT_mEc_Afc_AthetaC_psT_fDeltaH_RUT_a%% limit cycle parameters%% units are kmol, min, kJ, K, m^3%% 9/10/98, jbr%k_m=0.004;T_m=298;E=15000;c_Af=2;C_p=4;rho=1000;C_ps=C_p*rho;T_f=298;T_a=T_f;DeltaH_R=-2.2e5;U=340;theta=35;T_set=321.53;c_set=0.48995;T_fs=T_f;Kc=0;gamma=E/T_f;B=-DeltaH_R*c_Af*gamma/(C_ps*T_f);beta=U/C_ps*theta;Da=k_m*exp(-E*(1/T_f-1/T_m))*theta;x2c=(T_a-T_f)/T_f*gamma;%% building the lower branch;% use (theta,T) as dependent and% c_A as independent%x0=[1;T_f];nc_As=200;%tmp_table(1,:) = [0 T_f 0 -Inf -Inf 0];c_Avect=linspace(0.9999*c_Af,.005*c_Af,nc_As);fori=1:nc_Asc_A=c_Avect(i);opts=optimset('MaxFunEvals',2000*numel(x0),...'MaxIter',500*numel(x0));[x,fval,info]=fsolve('st_st_cA',x0,opts);theta=x(1);T=x(2);conv=(c_Af-c_A)/c_Af;%% compute the eigenvalues of the Jacobian%k=k_m*exp(-E*(1/T-1/T_m));Jac=[-1/theta-k,-k*c_A*E/(T*T);-k*DeltaH_R/C_ps,-U/C_ps-1/theta-k*c_A*DeltaH_R/C_ps*E/(T*T)];% lambda = sort(real(eig(Jac)))';lambda=eig(Jac);lamrow=[real(lambda(1))imag(lambda(1))real(lambda(2))...imag(lambda(2))];tmp_table(i,:)=[theta,T,conv,lamrow,info];x0=x;endtable=[tmp_table];%% turning the corner on the upper branch;% switch to (c_A,T) as dependent and% theta as independent% load x with current solution%x0=[c_A;T];cortheta=theta;nthetas=100;theta_vect=logspace(log10(cortheta),log10(500),nthetas);cleartmp_table;fori=1:nthetastheta=theta_vect(i);opts=optimset('MaxFunEvals',2000*numel(x0),...'MaxIter',500*numel(x0));[x,fval,info]=fsolve('st_st_theta',x0,opts);c_A=x(1);T=x(2);conv=(c_Af-c_A)/c_Af;k=k_m*exp(-E*(1/T-1/T_m));%% compute the eigenvalues of the Jacobian%Jac=[-1/theta-k,-k*c_A*E/(T*T);-k*DeltaH_R/C_ps,-U/C_ps-1/theta-k*c_A*DeltaH_R/C_ps*E/(T*T)];% lambda = sort(real(eig(Jac)))';lambda=eig(Jac);lamrow=[real(lambda(1))imag(lambda(1))real(lambda(2))...imag(lambda(2))];tmp_table(i,:)=[theta,T,conv,lamrow,info];x0=x;endtable=[table;tmp_table];save-asciist_st_osc.dattable;if(~strcmp(getenv('OMIT_PLOTS'),'true'))% PLOTTINGplot(table(:,1),table(:,3));axis([-5,100,0,1]);figure()semilogx(table(:,1),table(:,3));axis([.001,1000,0,1]);figure()plot(table(:,1),table(:,2));axis([-5,100,280,420]);figure()semilogx(table(:,1),table(:,2));axis([.001,1000,280,420]);% TITLEend% PLOTTING