% 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/16/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];npts1=200;xvect=logspace(log10(1e-4),log10(0.75),npts1);c_Avect=(1-xvect)*c_Af;cleartmp_table;fori=1:npts1c_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);a=real(lambda(1));b=imag(lambda(1));c=real(lambda(2));d=imag(lambda(2));if(a>=c)lamrow=[abcd];elselamrow=[cdab];endtmp_table(i,:)=[theta,T,conv,lamrow,info];x0=x;endtable=[tmp_table];%% stuff in some points near the extinction point to resolve% rapid change in the eigenvalues%npts2=50;xnow=conv;xvect=linspace(xnow,.95,npts2);c_Avect=(1-xvect)*c_Af;cleartmp_table;fori=1:npts2c_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;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=eig(Jac);a=real(lambda(1));b=imag(lambda(1));c=real(lambda(2));d=imag(lambda(2));if(a>=c)lamrow=[abcd];elselamrow=[cdab];endtmp_table(i,:)=[theta,T,conv,lamrow,info];x0=x;endtable=[table;tmp_table];%% after 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;npts3=200;theta_vect=logspace(log10(cortheta),log10(20.7),npts3);cleartmp_table;fori=1:npts3theta=theta_vect(i);[x,fval,info]=fsolve('st_st_theta',x0);c_A=x(1);T=x(2);conv=(c_Af-c_A)/c_Af;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=eig(Jac);a=real(lambda(1));b=imag(lambda(1));c=real(lambda(2));d=imag(lambda(2));if(a>=c)lamrow=[abcd];elselamrow=[cdab];endtmp_table(i,:)=[theta,T,conv,lamrow,info];x0=x;endtable=[table;tmp_table];%% jam a bunch of points in interval theta=[20.7,20.8] to pick up pair% of conjugate eigenvalues branching from real axis theta as% independent load x with current solution%x0=[c_A;T];cortheta=theta;npts4=100;theta_vect=linspace(cortheta,20.8,npts4);cleartmp_table;fori=1:npts4theta=theta_vect(i);[x,fval,info]=fsolve('st_st_theta',x0);c_A=x(1);T=x(2);conv=(c_Af-c_A)/c_Af;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=eig(Jac);a=real(lambda(1));b=imag(lambda(1));c=real(lambda(2));d=imag(lambda(2));if(a>=c)lamrow=[abcd];elselamrow=[cdab];endtmp_table(i,:)=[theta,T,conv,lamrow,info];x0=x;endtable=[table;tmp_table];%% finish out to large theta theta as independent load x with current% solution%x0=[c_A;T];cortheta=theta;npts5=200;theta_vect=logspace(log10(cortheta),log10(500),npts5);cleartmp_table;fori=1:npts5theta=theta_vect(i);[x,fval,info]=fsolve('st_st_theta',x0);c_A=x(1);T=x(2);conv=(c_Af-c_A)/c_Af;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=eig(Jac);a=real(lambda(1));b=imag(lambda(1));c=real(lambda(2));d=imag(lambda(2));if(a>=c)lamrow=[abcd];elselamrow=[cdab];endtmp_table(i,:)=[theta,T,conv,lamrow,info];x0=x;endtable=[table;tmp_table];%% for eigenvalue plot for s-part of curve, save rows 1:npts1+npts2%% for eigenvalue plot for limit cycle part of curve, save rows% 441:npts1+npts2+npts3+npts4+npts5%% Obviously the above row numbers change depending on how the points% were placed, so this part has to be redone if the problem parameters% change. I didn't see a better solution without going to a full% continuation approach, which seems like overkill for one or two% problems of this type. jbr 9/16/98rowend=npts1+npts2+npts3+npts4+npts5;tmp=table(441:rowend,:);save-asciist_st_osc_eigo.dattmp;if(~strcmp(getenv('OMIT_PLOTS'),'true'))% PLOTTINGplot(tmp(:,4),tmp(:,5),tmp(:,6),tmp(:,7));axis([-1,0.2,-.3,.3]);% TITLEend% PLOTTING