function test33
clear all;clc
L = 4;                  % ·´Ó¦¹Ü³¤, m
CA0=4.136;              % 2M1BĦ¶ûÁ÷ÂÊ, moles/hr
CB0=37.22  ;           % 2M2BĦ¶ûÁ÷ÂÊ, moles/hr
CC0 = 41.36;                  % MeOHĦ¶ûÁ÷ÂÊ, moles/hr
Cp = 131.33;                  % ±ÈÈÈ, kJ/kmol K
U = 0.105;                         % ´«ÈÈϵÊý, kJ/m2 hr K
dt = 0.025;                    % ¹Ü¾¶, m  
TJ = 333.15;                      % ÀäȴζÈ, K
T0 = 343.15;  
[z, y] = ode45(@Equations,[0 L],[CA0 CB0 CC0 T0]);
[z y(:,1) y(:,2) y(:,3) ]
CA = y(:, 1);
CB = y(:, 2);
CC = y(:, 3);
xA = (CA0-CA)./CA0;                        % AµÄת»¯ÂÊ
xB = (CB0-CB)./CB0;                    % BµÄת»¯ÂÊ
xB = [0; xB];
xC = (CC0-CC)./CC0;                    % CµÄת»¯ÂÊ
xC = [0; xC];

% ͼÐÎÊä³ö
figure(1)
plot(z, y(:, 4))        % ζȷֲ¼
xlabel('z')
ylabel('T (K)')
figure(2)
plot(z, xA, 'r-')       % ת»¯ÂÊ·Ö²¼
xlabel('z')
ylabel('x_A')
figure(3)
plot(z, CA, 'ro-', z, CB, 'm>-', z, CC, 'b*-')    % Ũ¶È·Ö²¼
xlabel('z')
ylabel('C_A, C_B, C_C')
legend('C_A', 'C_B', 'C_C')



function dydz = Equations(z, y)

rhoB = 770;             % ´ß»¯¼Á¶Ñ»ýÃܶÈ, kg/m3
rhog =691.82   ;         % ÒºÌå»ìºÏÎïÃܶÈ, kg/m3
CA0=4.136;              % 2M1BĦ¶ûÁ÷ÂÊ, moles/hr
CB0=37.22  ;           % 2M2BĦ¶ûÁ÷ÂÊ, moles/hr
CC0 = 41.36;                  % MeOHĦ¶ûÁ÷ÂÊ, moles/hr
Cp = 131.33;                  % ±ÈÈÈ, kJ/kmol K
U = 0.105;                         % ´«ÈÈϵÊý, kJ/m2 hr K
dt = 0.025;                    % ¹Ü¾¶, m  
TJ = 333.15;                      % ÀäȴζÈ, K
T0 = 343.15;                 % ÎïÁϽø¿ÚζÈ(³õʼζÈ), K
H1 = 33600;                  % ·´Ó¦2M1B¡úTµÄ·´Ó¦ÈÈ, kJ/kmol
H2 = 26800;             % ·´Ó¦2M2B¡úTµÄ·´Ó¦ÈÈ, kJ/kmol
us = 3.11   ;            % ¿ÕËÙ£¬h-1

% »î»¯ÄÜ, kJ/kmol
E1 =89830  ;           
E2 =120010 ;                                       

R = 8.314;                      % ÀíÏëÆøÌå³£Êý, kJ/kmol K

Ac = pi*(dt/2)^2;       % ·´Ó¦¹ÜµÄºá½ØÃæ»ý, m2



CA = y(1);                                    % A---2M1B
CB = y(2);                                    % B---2M2B
CC = y(3);                                    % C---MeOH
T = y(4);
% Ëٶȳ£Êý, kmol/kg catalyst hr
k1 = exp(-E1/(R*T) + 35.35);
k2 = exp(-E2/(R*T) + 42.47);
% ·´Ó¦Æ½ºâ³£Êý
K1= exp(-8.3881+4041.2/T);
K2= exp(-8.2473+3225.3/T);
% ·´Ó¦ËÙ¶È, kmol/kg catalyst hr   
rA = -k1*CA*CC + k1/K1*(CC0-CC)  ;       % 2M1BµÄÏûºÄËÙÂÊ
rB = -k2*CB*CC + k2/K2*(CC0-CC);        % 2M2BµÄÏûºÄËÙÂÊ
rC = -k1*CA*CC + k1/K1*(CC0-CC)-k2*CB*CC + k2/K2*(CC0-CC);   % ¼×´¼µÄ×Ü·´Ó¦ËÙ¶È

                         % Ä£ÐÍ·½³Ì×é


% ·´Ó¦ËÙ¶È


% ÎïÁÏƽºâ
dCAdz = rhoB*rA/us;
dCBdz = rhoB*rB/us;
dCCdz = rhoB*rC/us;

% ÈÈÁ¿ºâËã
dTdz = ( rhoB*((-H1*k1*y(1)*y(3))-(-H2*k2*y(2)*y(3)))-4*U*(T-TJ)/dt )/(us*rhog*Cp);
dydz = [dCAdz; dCBdz; dCCdz; dTdz];

2Â¥: Originally posted by ÔÂÖ»À¶ at 2014-02-17 21:02:13
ͼÐÎÈ總ͼ1¡£


function test33
clear all;clc
L = 4;                  % ·´Ó¦¹Ü³¤, m
CA0=4.136;              % 2M1BĦ¶ûÁ÷ÂÊ, moles/hr
CB0=37.22  ;           % 2M2BĦ¶ûÁ÷ÂÊ, moles/hr
CC0 = 41.3 ...

3Â¥: Originally posted by liqiang0915 at 2014-02-18 10:30:35
ËãµÄʱºò»¹ÊdzöÏÖ´íÎ󣬷³ÀÍÄúÔÙ¿´¿´¡£
??? Error using ==> feval
Undefined function or method 'Equations' for input arguments of type 'double'.

Error in ==> odearguments at 110
f0 = feval( ...