该楼层疑似违规已被系统折叠 隐藏此楼查看此楼
function test4
clc %清屏%
clear %消除内存变量%
[t,z]=ode45(@fd,[0:0.1:6],[0;753]); %用Matlab自带的库函数调用自己编写的方程组函数,求解微分方程组(1)(2)%
z1=z(:,1) %将z矩阵中的全部第一列数据赋给转化率z1%
z2=z(:,2) %将z矩阵中的全部第二列数据赋给床层温度z2%
figure(1),plot(t,z1,'r.-'), xlabel('L(m)'), ylabel('X_A') %绘制出转化率与催化剂床层高度之间的关系图, 其中l表示催化剂床层高度,z(;,1)表示转化率XA%
figure(2),plot(t,z2,'g.-'), xlabel('L(m)'), ylabel('T(K)')%绘制出床层温度与催化剂床层高度之间的关系图,其中z(;,2)表示床层温度T%
function dz=fd(t,z) %自己编写的微分方程组函数(1)(2)%
r=8.48*10e9*exp(-18410./z(2))*(12*(1-z(1))./(100+6.*z(1))); %输入化学反应速率已知条件表达式,单位kmol/(kg·cat·h)%
H=81500; %输入反应热参数,单位:KJ/kmol %
VA0=7066; %输入来气的体积流量,单位 m3/h %
VB0=16529; %输入稀释空气的体积流量,单位 m3/h %
rho1=680; %输入催化剂床层堆积密度,单位 kg/m3 %
rho2=1.301; %输入气体的平均密度,单位 kg/m3 %
cp1=10; %输入催化剂的比热,单位 KJ/(kg·K)%
cp2=1; %输入混合气体平均比热容,单位 KJ/(kg·K)%
D=1; %输入反应管内径,单位 m %
FA0=VA0*0.384/22.4; %计算输入的摩尔流量,单位 kmol/h %
u0=(VA0+VB0)/(3.14*(D/2)^2); %计算空塔气速,单位 m/h %
dz(1)=rho1*3.14*(D/2)^2*r./FA0; %输入方程(1)式的物料衡算已知条件式 %
dz(2)=rho1*H.*r./(u0*(rho1*cp1+rho2*cp2)); %输入方程(2)式的热量衡算已知条件式%
dz=[dz(1);dz(2)];