function IsothermTR
clear all;clc
global keq1 keq2 k1 k2 P
stoptime=1;dt=0.01;
T=573.9;keq1=0.312;keq2=0.480;
k1=1.496e7*exp(-8.44e3/T);k2=8.67e6*exp(-8.44e3/T);
P=1;F=128.2;
CA0=1/22.4;
CA0=CA0/453.6/(1e-3*3.2808^3);
v0=FA0/CA0;
x0=[0 0];t=[0,0.005,0.01,0.02,0.04,0.06,0.08,0.10,0.12,0.14,0.18,0.22,0.26,0.30,0.40];
[t,x]=ode45(@Equations,t,x0);
tau=CA0*t;SV=1./tau;
xBt=x(:,1)+x(:,2);
t42p=spline(xBt,t,0.42);
VR=t42p*FA0;
PB=1-x(:,1)-x(:,2);PD=x(:,1)/2-x(:,2);PT=x(:,2);
tmax=t(find(PD==max(PD)));
tau_max=tau(find(PD==max(PD)));
SV_max=SV(find(PD==max(PD)));
%输出结果
disp('result:')
fprintf('\n 所需的反应器体')
fprintf('t x1 x2,xBt PB PD PT')
disp([t x(:,1:2) xBt PB PD PT]),plot(t,xBt),xlabel('V/F_A_0,ft^3.hr/lbmol'),ylabel('x_B_t')
figure,plot(t,PB,'b--',t,PD,'m-',t,PT,'r-'),xlabel('V/F_A_0,ft^3.hr/lbmol'),ylabel('P,atm')
legend('P_B','P_D','P_T')
%----------------------------------------------
function dxdt=Equations(t,x)
global keq1 keq2 k1 k2 P
PB=P*(1-x(1)-x(2));PD=P*(x(1)/2-x(2));PH=P*(x(1)/2+x(2));PT=P*x(2);
r1=k1*(PB*PB-PD*PH/keq1);r2=k2*(PB*PD-PH*PT/keq2);
%---------------------------------------------
dx1dt=r1;dx2dt=r2;
dxdt=[dx1dt;dx2dt];