| 查看: 474 | 回复: 1 | ||
[求助]
程序调试
|
|
function PDEs2DS_Collocation % 用对称正交配置法求解固定床反应器二维拟均相稳态模型（二维稳态PDE方程组） % (只在r方向无因次化) % % Author: HUANG Huajiang % Copyright 2003 UNILAB Research Center, % East China University of Science and Technology, Shanghai, PRC % $Revision: 1.0 $ $Date: 2003/03/07 $ 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 CB0=4.314 ; % 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 H2 = -26800; % 反应2M2B→T的反应热, kJ/kmol us = 15.35 ; % 流速 m/h Dr=0.0001055 % m2/h G= 10337.76 %质量流速 kg/(m2h) % 活化能, kJ/kmol E =72610 ; R = 8.314; % 理想气体常数, kJ/kmol K Ac = pi*a^2; % 反应管的横截面积, m2 % Equation coefficient（方程的系数） L = 4; %反应管长度，m N = 3; % 内配置点个数 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'] % 求沿管长的平均转化率xa(i) for i=1:length(z) Cn = C(i, ;Ca(i) = quadl(@func,0,1,[],[],Rn,Cn)/(1^2/2) end % 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),C); 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 dydC = [dTdR dCdR]'; % ------------------------------------------------------------------ function rC = ReactionRate(T,C) % 计算反应速度 k = exp(-72610/(8.314*T) + 22.96); % 反应平衡常数 K= exp(-8.2473+3225.3/T); rC = -k.*C.^2 + k./K.*(4.314-C); % 2M2B的消耗速率 % ------------------------------------------------------------------ function y = func(R,Rn,Cn) C = spline(Rn,Cn,R); y = R.*x; |
回复此楼» 猜你喜欢
青椒八年已不青，大家都被折磨成啥样了？
已经有5人回复
-
青年基金C终止
已经有3人回复
-
26申博求博导推荐-遥感图像处理方向
已经有4人回复
-
限项规定
已经有7人回复
-
救命帖
已经有8人回复
-
招博士
已经有4人回复
-
西南交通大学国家级人才团队2026年博士研究生招生（考核制）—机械、材料、力学方向
已经有3人回复
-
英文综述是否需要润色及查重
已经有5人回复
-
为什么nbs上溴 没有产物点出现呢
已经有9人回复
-
最失望的一年
已经有18人回复
2楼2014-05-20 10:29:16