高等反应工程习题7有待改进.docx

1、高等反应工程习题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 求解过程计算框图求解问题：给

2、定进口温度(371)和出口一氧化碳含量(0.0212)，求催化床高度。2.1 逐步扫描法求催化床高度2.2 二分法求催化床高度3 编写计算程序，求出催化床高度运行环境：Matlab R2011a3.1 六个函数文件：函数文件(1)(5)在两个主程序中都有用到，函数文件(6)在二分法的主程序中用到。(1) rw.m% -Rg is the universal gas constant8.314J/K/molfunction rw1 = rw( Tb,yCO )pCO=yCO.*3.05; pCO2=(0.1298-yCO).*3.05; pH2=(0.4345-yCO).*3.05;pH2O=(

3、yCO+0.2925).*3.05;rw1=kT(Tb).*pCO.*pCO2.-0.5.*(1-pCO2.*pH2./(Kp(Tb).*pCO.*pH2O);%Calculate reaction rate,unit:kmol/kg/h function kT1 = kT(Tb) Rg=8.314; kT1=3.08.*10.6.*0.101325.-0.5.*exp(-104600./(Rg.*Tb);%Calculate kT end function Kp1=Kp(Tb)Kp1=exp(2.3026.*(2185./Tb-0.1102./2.3026.*log(Tb)+0.6218.*

4、10.-3.*Tb-1.0604.*10.-7.*Tb.2-2.218);%Calculate equilibrium constant Kp endend(2) HrCO.mfunction Hr = HrCO( Tb )Hr= -(10000+0.219.*Tb-2.845.*10.-3.*Tb.2+0.9703.*10.-6.*Tb.3).*4.184;%Calculate the heat of reaction end(3) cpM.mfunction cp = cpM( yCO,Tb )cp=yCO.*cpCO(Tb)+(yCO+0.2925).*cpH2O(Tb)+(0.1298

5、-yCO).*cpCO2(Tb)+(0.4345-yCO).*cpH2(Tb)+0.1432.*cpN2(Tb);%Calculate mixed gas isobaric molar heat capacity function cp1CO = 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

6、).*8.314; end function cp2H2O =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 function cp3CO2= cpCO2(Tb) pCO2= (0.1298-yCO).*3.05;cp3CO2=(3.18266+0.53754.*(Tb

7、./100)-0.020125.*(Tb./100).2+0.018520.*(pCO2./0.101325)-0.22009.*10.-2.*(pCO2./0.101325).*(Tb/100).*8.314; end function cp4H2= 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.18

8、4; end function cp5N2= 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; endend4) dyCOdl.m% -DR is the fixed bed diameter,unit:m% -pb is the catalyst bulk den

9、sity,unit:kg/m3% -COR is the catalyst correction coefficienthere COR is 0.380.% -NT0 is nlet gas flow,unit:kmol/h% -pi is circumference ratiofunction 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);%Calculate dyCO/dlend 5) dTbdl.mfunctio

10、n dTbdl1 = dTbdl( yCO,Tb )dTbdl1=HrCO(Tb)./cpM(yCO,Tb).*dyCOdl(yCO,Tb);%Calculate dTb/dlend6) RKS4yCO.mfunction e = RKS4yCO( L) yCO=;Tb=; yCO(1)=0.0810; Tb(1)=644.15; h=L./20; for i=1:20 k1=feval(dyCOdl,yCO(i),Tb(i); t1=feval(dTbdl,yCO(i),Tb(i);%Using Runge-Kuttas method to calculate yCO and Tb k2=f

11、eval(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

12、)=yCO(i)+h.*(k1+2.*k2+2.*k3+k4)./6; %Calculate exit yCO of the given lenth Tb(i+1)=Tb(i)+h.*(t1+2.*t2+2.*t3+t4)./6; end e=yCO(i+1)-0.0212;end3.2 主程序与运行结果3.2.1 逐步扫描法主程序：calculateyCOTb.m% -yCO is the carbon monoxide molar fraction% -Tb is the bed temperature,unit:K% -DeltHrCO is the heat of reaction,u

13、nit:J/mol% -cpMix is the mixed gas isobaric molar heat capacity, unit:J/K/mol% -L is the height of the bed,unit:myCO=;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;while yCO(i)-0.02120.00001 k1=feval(dyCOdl,yCO(i),Tb(i); t1=f

14、eval(dTbdl,yCO(i),Tb(i);%Using Runge-Kuttas method to calculate yCO and Tb 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,(yC

15、O(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;endR=L,(Tb-273.15),Tb,yCO,DeltHrCO,cpMix;disp(-)format short g dis

16、p( 床高 床层温度 床层温度 一氧化碳分率 反应热 混合气体等压摩尔热容)disp( l/m tb/ tb/K yCO DeltHr/(J/mol) CpMix/(J/mol/K)disp(-)disp(R)disp(-)主程序运行结果： 可以发现，床层温度为4.32m时，出口一氧化碳摩尔分率达到0.0212，满足给定要求。3.2.2 二分法通过试算，发现当床层高度为2m时，当床层高度为5m时，因此，选取床层高度的计算范围为2,5。主程序：dichotomycalculateL.ma=2;b=5;while abs(b-a)0.0001 x=(b+a)./2; if RKS4yCO(x)=0

17、 break elseif sign( RKS4yCO(x)=sign(RKS4yCO(b) b=x; else a=x; endend yCO=;Tb=; yCO(1)=0.0810; Tb(1)=644.15; h=x./20;DR=4; for i=1:20 k1=feval(dyCOdl,yCO(i),Tb(i); t1=feval(dTbdl,yCO(i),Tb(i);%Using Runge-Kuttas method to calculate yCO and Tb k2=feval(dyCOdl,(yCO(i)+h./2.*k1),(Tb(i)+h./2.*t1); t2=fev

18、al(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.*(

19、t1+2.*t2+2.*t3+t4)./6; DeltHrCO(i+1)=HrCO(Tb(i+1); cpMix(i+1)=cpM(yCO(i),Tb(i); endR=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，满足给定要求。二分法与逐步扫描法的

20、运算结果相近。4 给定出口一氧化碳含量(0.0212)，确定最佳进口温度（选做）试算时，发现进口温度超过一定数值时，龙格库塔法求解yCO和Tb过程中的k1和t1的值极小，导致计算陷入死循环。改变进口温度并经过多次试算后，确定进口温度计算范围的终点取660.15K，进口温度的始点这里取638.15K。由于催化剂体积对进口温度的一阶微分方程不容易获得，本题的求解思路是通过改变进口温度得到相应的催化剂体积，通过对运算结果列表的观察确定最佳进口温度。 3.1的中函数文件(1)(5)在下面的主程序中有用到。主程序：calculatebestTb.myCO=;Tb=;L=;V=;Lenth=;Tin=;T

21、out=;yCOin=;yCOout=;yCO(1)=0.0810; L(1)=0;h=0.02;j=1;for T=638.15:1:660.15 Tb(1)=T;i=1;DR=4;while yCO(i)-0.02120.00001 k1=feval(dyCOdl,yCO(i),Tb(i); t1=feval(dTbdl,yCO(i),Tb(i);%Using Runge-Kuttas method to calculate yCO and Tb k2=feval(dyCOdl,(yCO(i)+h./2.*k1),(Tb(i)+h./2.*t1); t2=feval(dTbdl,(yCO(

22、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

23、+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;endformat short gR=(Tin-273.15),(Tout-273.15),yCOin,yCOout,Lenth,V;disp(-)disp( 进口温度 出口温度 进口一氧化碳 出口一氧化碳 床层高度 催化剂体积(COR为0.380)disp( / / 摩尔分率 摩尔分率 /m /m3 )disp(-)disp(R)disp(-)主程序的运算结果： 可以发现，在进口温度为365387的范围内，催化剂体积随进口温度的升高而下降。当催化剂体积当进口温度为383时，催化剂的体积最小，为46.998m3。