|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
function PDEs2DS_324
clear all
clc
global r Fco Fco2 F H1 H2 Cpm N u L kg1a pb pa kg2a
% Parameters
r = 1; % W/(m K)
Fco = 14; % kg/(m^2 h)
Fco2 = 4; % kJ/(kg K)
F = 200; % J/mol
H1=206100; % kg/m^3
H2=164700; % radius of reactor, m
Cpm =31.02577; % u0*c0 obtained from u0*c0*A=0.069 kmol/h
N=11; % CP', J/(kg K)
u=1.4385; % kg/s
L=3;
kg1a=6.76;
pb=1175;
pa=2;
kg2a=5.83;
y0 = [543.3 574.5203462 607.6242169 640.1415255 671.1218255 700.2984569 727.6492903 753.2367257 777.1515236 799.4921566 820.3568002 0 0.055093867 0.113920607 0.171932185 0.227264584 0.279360997 0.328151221 0.373737811 0.416283021 0.455966554 0.492968961 0 0.052761405 0.106914317 0.15911271 0.208568428 0.255208536 0.299132063 0.340477585 0.379389327 0.416007726 0.450466717]; % y=[T1 T2 T3 T4 x1 x2 x3]
tspan=0:1:3000;
[t,y]=ode45(@Euqations,tspan,y0);
T = y(:,1:N);
x1 = y(:,N+1:2*N);
x2 = y(:,2*N+1:3*N);
% ------------------------------------------------------------------
function dydx = Euqations(x1,x2,y)
global r Fco Fco2 F H1 H2 Cpm N u L
detaz=L/(N-1);
T = y(1:N);
x1 = y(N+1:2*N);
x2 = y(2*N+1:3*N);
rco= ReactionRate1(T(1:N),x1);
rco2=ReactionRate2(T(1:N),x2);
dTdt(1) = u*((T(1)-T(2))/detaz+(pi*r*r/F)*(rco(1)*H1+rco2(1)*H2)/Cpm);
dx1dt(1) =u*((x1(1)-x1(2))/detaz+(pi*r*r/Fco)*rco(1));
dx2dt(1) =u*((x2(1)-x2(2))/detaz+(pi*r*r/Fco2)*rco2(1));
dTdt(N) = u*((T(N-1)-T(N))/detaz+(pi*r*r/F)*(rco(N)*H1+rco2(N)*H2)/Cpm);
dx1dt(N) = u*((x1(N-1)-x1(N))/detaz+(pi*r*r/Fco)*rco(N));
dx2dt(N) = u*((x1(N-1)-x1(N))/detaz+(pi*r*r/Fco2)*rco(N));
for i = 2:N-1
dTdt(i) = u*((T(i-1)-T(i+1))/2/detaz+(pi*r*r/F)*(rco(i)*H1+rco2(i)*H2)/Cpm);
dx1dt(i) =u*((x1(i-1)-x1(i+1))/2/detaz+(pi*r*r/Fco)*rco(i));
dx2dt(i) =u*((x1(i-1)-x1(i+1))/2/detaz+(pi*r*r/Fco)*rco(i));
end
dydx = [dTdR dx1dt dx2dt]';
% ------------------------------------------------------------------
function f1 = ReactionRate1(T,x,T0) % 计算反应速度
global kg1a pb pa
CA0=pa*100000/8.314/T0*0.07;
f1=(kg1a*CA0*(1-x)*pb*533488.642.*exp(-76856.4/8.314./T).*(CA0*(1-x)).^0.673833)/20./(kg1a*CA0*(1-x)+pb*533529*exp(-76857/8.314./T).*(CA0*(1-x)).^0.674);
function f2 = ReactionRate2(T,x,T0) % 计算反应速度
global kg2a pb pa
CB0=pa*100000/8.314/T0*0.02;
f2=(kg2a*CB0*(1-x)*pb*296452638.*exp(-105858/8.314./T).*(CB0*(1-x)).^0.008)/20./(kg2a*CB0*(1-x)+pb*296452638*exp(-105858/8.314./T).*(CB0*(1-x)).^0.008);
% ------------------------------------------------------------------
|
|