Input argument "y" is undefined.请各位帮下忙
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 = ; % y= tspan=0:1:3000; =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 = '; % ------------------------------------------------------------------ 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); % ------------------------------------------------------------------ 1# lxx244lxx程序在下面，上面太乱了
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 = ;% y=
tspan=0:1:3000;
=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 = ';
% ------------------------------------------------------------------
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);
% ------------------------------------------------------------------ 你的方程定义的有问题，x1和x2可不可以放在一个数组中？
页:
[1]