24小时热门版块排行榜

>论坛更新日志 (4153)
>文献求助 (435)
>导师招生 (380)
>虫友互识 (373)
>休闲灌水 (255)
>考博 (171)
>硕博家园 (125)
>招聘信息布告栏 (120)
>博后之家 (98)
>论文投稿 (91)
>公派出国 (75)
>基金申请 (73)
>考研 (68)
>论文道贺祈福 (57)
>找工作 (54)
>教师之家 (41)

返回列表

liqiang0915

木虫 (初入文坛)

应助: 0 (幼儿园)
金币: 1718
帖子: 18
在线: 352.6小时
虫号: 1455044
注册: 2011-10-22
专业: 化学反应工程

[求助] Matlab求解固定床二维均相模型的代码，请Matlab大神帮忙看看已有1人参与

function PDEs2DS_Collocation
% 用对称正交配置法求解固定床反应器二维拟均相稳态模型（二维稳态PDE方程组）
% (只在r方向无因次化)

clear all
clc

global A B a  Rn T0 L N

% P3(x^2)对于圆柱对称（a=2）的配置常数Aji和Bji
A = [-3.359794  5.2924315 -3.1010284  1.1683909
   -1.3980385 -1.5627540  4.3197367 -1.3589422
   0.69721650 -3.6766754  -1.1267583  4.1062172
   -1.2266754 5.4010626 -19.174383  15];
B = [-15.881426 19.636380 -5.2811862  1.5262327
   11.151861    -34.497415  29.235709 -5.890155
   -3.5405872 34.512110 -99.621159  68.649637
   -33.869987 136.21969 -252.37970  150];

Rn = [0.29763730  0.63989598  0.88750181  1]; % 全部配置点

% Parameters
rho = 607.3;          % 催化剂堆积密度, kg/m3
rhog =673.47 ;       % 液体混合物密度, kg/m3
CA0=0.963;          % 2M1B摩尔流率, kmoles/m3 hr
CB0=3.351  ;          % 2M2B摩尔流率, kmoles/m3 hr
CC0 = 4.314;          % MeOH摩尔流率, kmoles/m3 hr
Cp = 2.828;                % 比热, kJ/kg K
Ramda = 364.32;             % 传热系数, kJ/m2 hr K
a= 0.0125;                   % 管半径, m  radius of reactor, m
TJ = 333.15;                   % 冷却温度, K
T0 = 333.15;                % 物料进口温度(初始温度), K
H1 = -33600;                % 反应2M1B→T的反应热,kJ/kmol
H2 = -26800;          % 反应2M2B→T的反应热, kJ/kmol
us = 15.35 ;          % 流速 m/h
Dr=0.0001055             % m2/h
G= 10337.76             %质量流速 kg/(m2h)

% 活化能, kJ/kmol
E1 =85090  ;
E2 =72610 ;

R = 8.314;             % 理想气体常数, kJ/kmol K

Ac = pi*a^2;    % 反应管的横截面积, m2

% Equation coefficient（方程的系数）
L = 4;             %反应管长度，m
N = 3;             % 内配置点个数

[z, q] = ode45(@Equations,[0 L],[CA0 CB0 CC0 T0]);
[z q(:,1) q(:,2) q(:,3) q(:,4) ]
CA = q(:, 1);
CB = q(:, 2);
CC = q(:, 3);
T = q(:,4);

y0 = [T0 T0 T0 333.15 0 0 0];  % y=[T1 T2 T3 T4 C1 C2 C3]
[z,y] = ode45(@Euqations,[0 L],y0);
z
T = y(:,1:N+1)
C = y(:,N+2:2*N+1);    % TAME浓度
for i = 1:length(z)
Cb(i) = - sum(A(N+1,1:N).* C(i,

/A(N+1,N+1));
end
C = [C Cb']

% Plot the results
surf(Rn*a,z,T)    % 反应管轴径向温度分布
xlabel('r (m)')
ylabel('z (m)')
zlabel('T (K)')
figure
plot(z,Ca)       % 平均转化率沿管长的分布图
xlabel('z (m)')
ylabel('C_a_v')
figure
surf(Rn*a,z,C)    % 轴径向平均转化率分布
xlabel('r (m)')
ylabel('z (m)')
zlabel('C')

% ------------------------------------------------------------------
function dydC = Euqations(C,y)
global A B a  Rn T0 L N
T = y(1:N+1);
C = y(N+2:2*N+1);
rC = ReactionRate(T(1:N),CA,CB,CC);
for i = 1:N
dTdR(i) =Ramda*L/Cp/G/a^2* sum( (B(i,

+A(i,

./Rn(i)) .* T' ) - rho*rC(i);
end

dTdR(N+1) = sum(A(N+1,

.* T')

Cb = - sum(A(N+1,1:N).* C'/A(N+1,N+1));
for i=1:N
dCdR(i) = Dr*L/u/a^2*( sum((B(i,1:N)+A(i,1:N)./Rn(i)).*C(1:N)') ...
      + (B(i,N+1)+A(i,N+1)./Rn(i)).*Cb ) + rho*rC(i);
end
dydx = [dTdR dCdR]';

% ------------------------------------------------------------------
function rC = ReactionRate(T,CA,CB,CC)    % 计算反应速度
k1 = exp(-E1/(R*T) + 30.26);
k2 = exp(-E2/(R*T) + 22.96);
% 反应平衡常数
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 = rho*rA/us;
dCBdz = rho*rB/us;
dCCdz = rho*rC/us;

回复此楼