高等反应工程习题7有待改进.docx
《高等反应工程习题7有待改进.docx》由会员分享，可在线阅读，更多相关《高等反应工程习题7有待改进.docx（15页珍藏版）》请在冰豆网上搜索。
高等反应工程习题7有待改进
高等反应工程习题7（例4-1）
本小组（第7组）的催化剂活性校正系数COR取0.380。
1单段固定催化床一维拟均相模型
反应速率方程：
（1-1）
热量衡算式：
（1-2）
对于轴向的反应微元，有
（1-3）
（1-4）
对于一氧化碳变换反应，有
（1-5）
此外，
（1-6）
由式（1-1）~式（1-6），经过整理可以得到
（1-7）
（1-8）
几个重要公式：
本征速率方程
（1-9）
式（1-9）中，
（1-10）
（1-10）
反应热
（1-11）
混合气体等压摩尔热容
（1-12）
（1-13）
（1-14）
（1-15）
（1-16）
（1-17）
2求解过程计算框图
求解问题：
给定进口温度（371℃）和出口一氧化碳含量（0.0212），求催化床高度。
2.1逐步扫描法求催化床高度
2.2二分法求催化床高度
3编写计算程序，求出催化床高度
运行环境：
MatlabR2011a
3.1六个函数文件：
函数文件
（1）~（5）在两个主程序中都有用到，函数文件（6）在二分法的主程序中用到。
（1）rw.m
%-Rgistheuniversalgasconstant£¬8.314J/K/mol
function[rw1]=rw（Tb,yCO）
pCO=yCO.*3.05;pCO2=（0.1298-yCO）.*3.05;pH2=（0.4345-yCO）.*3.05;pH2O=（yCO+0.2925）.*3.05;
rw1=kT（Tb）.*pCO.*pCO2.^-0.5.*（1-pCO2.*pH2./（Kp（Tb）.*pCO.*pH2O））;%Calculatereactionrate,unit:
kmol/kg/h
functionkT1=kT（Tb）
Rg=8.314;
kT1=3.08.*10.^6.*0.101325.^-0.5.*exp（-104600./（Rg.*Tb））;%CalculatekT
end
functionKp1=Kp（Tb）
Kp1=exp（2.3026.*（2185./Tb-0.1102./2.3026.*log（Tb）+0.6218.*10.^-3.*Tb-1.0604.*10.^-7.*Tb.^2-2.218））;%CalculateequilibriumconstantKp
end
end
（2）HrCO.m
function[Hr]=HrCO（Tb）
Hr=-（10000+0.219.*Tb-2.845.*10.^-3.*Tb.^2+0.9703.*10.^-6.*Tb.^3）.*4.184;%Calculatetheheatofreaction
end
（3）cpM.m
function[cp]=cpM（yCO,Tb）
cp=yCO.*cpCO（Tb）+（yCO+0.2925）.*cpH2O（Tb）+（0.1298-yCO）.*cpCO2（Tb）+（0.4345-yCO）.*cpH2（Tb）+0.1432.*cpN2（Tb）;%Calculatemixedgasisobaricmolarheatcapacity
functioncp1CO=cpCO（Tb）
pCO=yCO.*3.05;
cp1CO=（3.86771-0.23279.*（Tb/100）+0.046135.*（Tb/100）.^2-0.2168.*10.^-2.*（Tb/100）.^3+0.42112.*10.^-2.*（pCO./0.101325）-0.4694.*10^-3.*（pCO./0.101325）.*（Tb./100））.*8.314;
end
functioncp2H2O=cpH2O（Tb）
pH2O=（yCO+0.2925）.*3.05;
cp2H2O=（0.65765-0.049712.*（Tb/100）+0.5269.*10.^-3.*（Tb./100）.^-3+0.020739.*（pH2O./0.101325）-0.27123.*10^-2.*（pH2O./0.101325）.*（Tb./100））.*18.0513.*4.184;
end
functioncp3CO2=cpCO2（Tb）
pCO2=（0.1298-yCO）.*3.05;
cp3CO2=（3.18266+0.53754.*（Tb./100）-0.020125.*（Tb./100）.^2+0.018520.*（pCO2./0.101325）-0.22009.*10.^-2.*（pCO2./0.101325）.*（Tb/100））.*8.314;
end
functioncp4H2=cpH2（Tb）
pH2=（0.4345-yCO）.*3.05;
cp4H2=（6.8712+0.03135.*（Tb/100）+0.14138.*10.^-2.*（pH2./0.101325）-0.6.*10.^-6.*（pH2./0.101325）.^2+0.1603.*10.^-3.*（pH2./0.101325）.*（Tb/100））.*4.184;
end
functioncp5N2=cpN2（Tb）
pN2=0.1432.*3.05;
cp5N2=（4.23329-0.4145.*（Tb./100）+0.072309.*（Tb./100）.^2-0.34116.*10.^-2.*（Tb./100）.^3+0.57726.*10.^-2.*（pN2./0.101325）-0.7404.*10.^-3.*（pN2./0.101325）.*（Tb./100））.*8.314;
end
end
4）dyCOdl.m
%-DRisthefixedbeddiameter,unit:
m
%-pbisthecatalystbulkdensity,unit:
kg/m3
%-CORisthecatalystcorrectioncoefficient£¬hereCORis0.380.
%-NT0isnletgasflow,unit:
kmol/h
%-piiscircumferenceratio
function[dyCOdl1]=dyCOdl（yCO,Tb）
DR=4;pb=1500;COR=0.380;NT0=9707.4;
dyCOdl1=pi./4.*DR.^2.*pb.*COR./NT0.*（-rw（Tb,yCO））;%CalculatedyCO/dl
end
5）dTbdl.m
function[dTbdl1]=dTbdl（yCO,Tb）
dTbdl1=HrCO（Tb）./cpM（yCO,Tb）.*dyCOdl（yCO,Tb）;%CalculatedTb/dl
end
6）RKS4yCO.m
function[e]=RKS4yCO（L）
yCO=[];Tb=[];
yCO
（1）=0.0810;Tb
（1）=644.15;
h=L./20;
fori=1:
20
k1=feval（'dyCOdl',yCO（i）,Tb（i））;t1=feval（'dTbdl',yCO（i）,Tb（i））;%UsingRunge-Kutta'smethodtocalculateyCOandTb
k2=feval（'dyCOdl',（yCO（i）+h./2.*k1）,（Tb（i）+h./2.*t1））;t2=feval（'dTbdl',（yCO（i）+h./2.*k1）,（Tb（i）+h./2.*t1））;
k3=feval（'dyCOdl',（yCO（i）+h./2.*k2）,（Tb（i）+h./2.*t2））;t3=feval（'dTbdl',（yCO（i）+h./2.*k2）,（Tb（i）+h./2.*t2））;
k4=feval（'dyCOdl',（yCO（i）+h.*k3）,（Tb（i）+h.*t3））;t4=feval（'dTbdl',（yCO（i）+h.*k3）,（Tb（i）+h.*t3））;
yCO（i+1）=yCO（i）+h.*（k1+2.*k2+2.*k3+k4）./6;%CalculateexityCOofthegivenlenth
Tb（i+1）=Tb（i）+h.*（t1+2.*t2+2.*t3+t4）./6;
end
e=yCO（i+1）-0.0212;
end
3.2主程序与运行结果
3.2.1逐步扫描法
主程序：
calculateyCOTb.m
%-yCOisthecarbonmonoxidemolarfraction
%-Tbisthebedtemperature,unit:
K
%-DeltHrCOistheheatofreaction,unit:
J/mol
%-cpMixisthemixedgasisobaricmolarheatcapacity,unit:
J/K/mol
%-Listheheightofthebed,unit:
m
yCO=[];Tb=[];L=[];DeltHrCO=[];cpMix=[];
yCO
（1）=0.0810;Tb
（1）=644.15;L
（1）=0;h=0.02;
DeltHrCO
（1）=HrCO（Tb
（1））;cpMix
（1）=cpM（yCO
（1）,Tb
（1））;
i=1;
whileyCO（i）-0.0212>0.00001
k1=feval（'dyCOdl',yCO（i）,Tb（i））;t1=feval（'dTbdl',yCO（i）,Tb（i））;%UsingRunge-Kutta'smethodtocalculateyCOandTb
k2=feval（'dyCOdl',（yCO（i）+h./2.*k1）,（Tb（i）+h./2.*t1））;t2=feval（'dTbdl',（yCO（i）+h./2.*k1）,（Tb（i）+h./2.*t1））;
k3=feval（'dyCOdl',（yCO（i）+h./2.*k2）,（Tb（i）+h./2.*t2））;t3=feval（'dTbdl',（yCO（i）+h./2.*k2）,（Tb（i）+h./2.*t2））;
k4=feval（'dyCOdl',（yCO（i）+h.*k3）,（Tb（i）+h.*t3））;t4=feval（'dTbdl',（yCO（i）+h.*k3）,（Tb（i）+h.*t3））;
yCO（i+1）=yCO（i）+h.*（k1+2.*k2+2.*k3+k4）./6;Tb（i+1）=Tb（i）+h.*（t1+2.*t2+2.*t3+t4）./6;
L（i+1）=L（i）+h;
DeltHrCO（i+1）=HrCO（Tb（i+1））;
cpMix（i+1）=cpM（yCO（i）,Tb（i））;
i=i+1;
end
R=[L',（Tb'-273.15）,Tb',yCO',DeltHrCO',cpMix'];
disp（'------------------------------------------------------------------------------------------'）
formatshortg
disp（'床高´床层温度床层温度一氧化碳分率反应热混合气体等压摩尔热容'）
disp（'l/mtb/℃tb/KyCODeltHr/（J/mol）CpMix/（J/mol/K）'）
disp（'------------------------------------------------------------------------------------------'）
disp（R）
disp（'------------------------------------------------------------------------------------------'）
主程序运行结果：
…
可以发现，床层温度为4.32m时，出口一氧化碳摩尔分率达到0.0212，满足给定要求。
3.2.2二分法
通过试算，发现当床层高度为2m时，
，当床层高度为5m时，
，因此，选取床层高度的计算范围为[2,5]。
主程序：
dichotomycalculateL.m
a=2;b=5;
whileabs（b-a）>0.0001
x=（b+a）./2;
ifRKS4yCO（x）==0
break
elseifsign（RKS4yCO（x））==sign（RKS4yCO（b））
b=x;
else
a=x;
end
end
yCO=[];Tb=[];
yCO
（1）=0.0810;Tb
（1）=644.15;
h=x./20;DR=4;
fori=1:
20
k1=feval（'dyCOdl',yCO（i）,Tb（i））;t1=feval（'dTbdl',yCO（i）,Tb（i））;%UsingRunge-Kutta'smethodtocalculateyCOandTb
k2=feval（'dyCOdl',（yCO（i）+h./2.*k1）,（Tb（i）+h./2.*t1））;t2=feval（'dTbdl',（yCO（i）+h./2.*k1）,（Tb（i）+h./2.*t1））;
k3=feval（'dyCOdl',（yCO（i）+h./2.*k2）,（Tb（i）+h./2.*t2））;t3=feval（'dTbdl',（yCO（i）+h./2.*k2）,（Tb（i）+h./2.*t2））;
k4=feval（'dyCOdl',（yCO（i）+h.*k3）,（Tb（i）+h.*t3））;t4=feval（'dTbdl',（yCO（i）+h.*k3）,（Tb（i）+h.*t3））;
yCO（i+1）=yCO（i）+h.*（k1+2.*k2+2.*k3+k4）./6;
Tb（i+1）=Tb（i）+h.*（t1+2.*t2+2.*t3+t4）./6;
DeltHrCO（i+1）=HrCO（Tb（i+1））;
cpMix（i+1）=cpM（yCO（i）,Tb（i））;
end
R=[371,Tb（i+1）,0.0810,yCO（i+1）,x,pi./4.*DR.^2.*x];
disp（'-------------------------------------------------------------------------------------------'）
disp（'进口温度出口温度进口一氧化碳出口一氧化碳床层高度催化剂体积（COR为0.380）'）
disp（'/℃/℃摩尔分率摩尔分率/m/m3'）
disp（'-------------------------------------------------------------------------------------------'）
disp（R）
disp（'-------------------------------------------------------------------------------------------'）
主程序运行结果：
可以发现，当床层高度为4.3206m时，出口一氧化碳摩尔分率为0.0212，满足给定要求。
二分法与逐步扫描法的运算结果相近。
4给定出口一氧化碳含量（0.0212），确定最佳进口温度（选做）
试算时，发现进口温度超过一定数值时，龙格库塔法求解yCO和Tb过程中的k1和t1的值极小，导致计算陷入死循环。
改变进口温度并经过多次试算后，确定进口温度计算范围的终点取660.15K，进口温度的始点这里取638.15K。
由于催化剂体积对进口温度的一阶微分方程不容易获得，本题的求解思路是通过改变进口温度得到相应的催化剂体积，通过对运算结果列表的观察确定最佳进口温度。
3.1的中函数文件
（1）~（5）在下面的主程序中有用到。
主程序：
calculatebestTb.m
yCO=[];Tb=[];L=[];
V=[];Lenth=[];Tin=[];Tout=[];yCOin=[];yCOout=[];
yCO
（1）=0.0810;L
（1）=0;h=0.02;
j=1;
forT=638.15:
1:
660.15
Tb
（1）=T;i=1;DR=4;
whileyCO（i）-0.0212>0.00001
k1=feval（'dyCOdl',yCO（i）,Tb（i））;t1=feval（'dTbdl',yCO（i）,Tb（i））;%UsingRunge-Kutta'smethodtocalculateyCOandTb
k2=feval（'dyCOdl',（yCO（i）+h./2.*k1）,（Tb（i）+h./2.*t1））;t2=feval（'dTbdl',（yCO（i）+h./2.*k1）,（Tb（i）+h./2.*t1））;
k3=feval（'dyCOdl',（yCO（i）+h./2.*k2）,（Tb（i）+h./2.*t2））;t3=feval（'dTbdl',（yCO（i）+h./2.*k2）,（Tb（i）+h./2.*t2））;
k4=feval（'dyCOdl',（yCO（i）+h.*k3）,（Tb（i）+h.*t3））;t4=feval（'dTbdl',（yCO（i）+h.*k3）,（Tb（i）+h.*t3））;
yCO（i+1）=yCO（i）+h.*（k1+2.*k2+2.*k3+k4）./6;Tb（i+1）=Tb（i）+h.*（t1+2.*t2+2.*t3+t4）./6;
L（i+1）=L（i）+h;
DeltHrCO（i+1）=HrCO（Tb（i+1））;
cpMix（i+1）=cpM（yCO（i）,Tb（i））;
i=i+1;
end
Tin（j）=T;Tout（j）=Tb（i）;
yCOin（j）=0.0810;yCOout（j）=yCO（i）;
Lenth（j）=L（i）;V（j）=pi./4.*DR.^2.*L（i）;
j=j+1;
end
formatshortg
R=[（Tin'-273.15）,（Tout'-273.15）,yCOin',yCOout',Lenth',V'];
disp（'------------------------------------------------------------------------------------------'）
disp（'进口温度出口温度进口一氧化碳出口一氧化碳床层高度催化剂体积（COR为0.380）'）
disp（'/℃/℃摩尔分率摩尔分率/m/m3'）
disp（'------------------------------------------------------------------------------------------'）
disp（R）
disp（'------------------------------------------------------------------------------------------'）
主程序的运算结果：
可以发现，在进口温度为365~387℃的范围内，催化剂体积随进口温度的升高而下降。
当催化剂体积当进口温度为383℃时，催化剂的体积最小，为46.998m3。
升级会员 