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
算的时候还是出现错误，烦劳您再看看。
??? Error using ==> feval
Undefined function or method 'Equations' for input arguments of type 'double'.

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