| 查看: 1143 | 回复: 1 | ||
[求助]
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; |
回复此楼» 猜你喜欢
瑞典林雪平大学博士后招聘｜PEC Water Splitting 方向
已经有130人回复
-
PRB投稿，Acknowledgment sent to author状态十天了
已经有3人回复
-
物理学I论文润色/翻译怎么收费?
已经有154人回复
-
PVA溶解
已经有2人回复
-
基金申请
已经有44人回复
-
CSC与新西兰维多利亚大学PhD奖学金项目
已经有0人回复
-
新西兰Robinson研究所 招聘CSC公派访问人员
已经有0人回复
-
帮我的英语口语老师找学生
已经有0人回复
-
什么时候开奖？
已经有13人回复
-
散金币祈福
已经有94人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- matlab动力学模型编程求助
已经有6人回复
- 求回归方程参数的最优化算法matlab代码
已经有17人回复
- 求大神用matlab给我处理一组数据，有人没有啊？
已经有15人回复
- 中心差分求解二维热传导的matlab程序
已经有12人回复
- 跪求永磁同步风力发电机组在ＭＡＴＬＡＢ搭建的仿真模型
已经有5人回复
- 各位大神，请问我的Matlab批处理程序怎么改才能用exce一次l输出所有数据
已经有3人回复
- 关于matlab进行傅里叶变换和逆变换的程序，求大神指导！
已经有9人回复
- 求基于DCT域的数字图像水印Matlab代码，谢谢各位大神！
已经有3人回复
- matlab解决重心法选址问题，运行不了啊，求助各位大神
已经有10人回复
- 请教MATLAB中的LMI求解问题，急！！！
已经有7人回复
- MATLAB微分方程参数拟合问题，求大神
已经有7人回复
- 请教大家一个反卷积的问题~~帮忙推导或者matlab编程计算~~
已经有5人回复
- 求助matlab一程序运行的问题，帮忙看看哪里不对
已经有4人回复
- 帮忙看看这个平面图用matlab怎么画
已经有19人回复
- matlab作业哪位大神麻烦给做一下,十万火急,周五之间要交
已经有3人回复
- 求助啊！求一程序，用matlab程序做，用蒙特卡洛方法模拟
已经有10人回复
- 【求助】请教虫友关于matlab差分方程的求解和画图程序
已经有4人回复
- 【求助】 求MAtlab求解程序！！！
已经有20人回复
lmg0608
铁杆木虫 (知名作家)
- 应助: 0 (幼儿园)
- 金币: 10621.2
- 散金: 30
- 红花: 1
- 帖子: 7439
- 在线: 469.3小时
- 虫号: 3647756
- 注册: 2015-01-16
- 专业: 传热传质学
2楼2017-07-26 09:48:25
28