1、二氧化硫轉化率最優化化研1313 王曉冰 2013200186SO21/2O2=SO3，四段絕熱反應器，級間間接換熱。1.基礎數據：混合物恒壓熱容Cp 0.2549 kcal/kgKH =23135 kcal/kmol催化劑堆密度b554 kg/m3進口SO2濃度8.0%，O2濃度9.0%，其余為氮氣處理量131 kmol SO2/hr，要求最終轉化率98 2.動力學方程：式中：3.基本要求：(1)在TX圖上，做出平衡線，至少4條等速率線；(2)以一維擬均相平推流模型為基礎，在催化劑用量最少的前提下，總的及各段的催化劑裝量；進出口溫度、轉化率；并在T-X圖上標出折線；4.討論：(1)要求的最終

2、轉化率從98變化到99對催化劑用量的影響；(2)YO2YSO221，SO2進口濃度在79之間變化，對催化劑裝量的影響。計算過程：編寫c+程序，如下：#include<math.h>#include<stdio.h>#define R 1.987#define h 0.00001double r(double x, double t, double xso2)double r,r1,r2,r3, B,keff,K, Kp,Pso2,Pso3,Po2; if(t>=693.15 && t<748.15)keff = 7.6915*pow(10,1

3、8)*exp(-76062 / (R*t); if(t>=748.15 && t<=873.15) keff = 1.5128*pow(10,7)*exp(-35992 / (R*t); K = 2.3*pow(10,-8)*exp(27200/(R*t); Kp=2.26203*pow(10,-5)*exp(11295.3 / t); Pso2=(xso2-xso2*x) / (1-xso2*x/2); Pso3=xso2*x / (1-xso2*x/2); Po2=(0.17-xso2-xso2*x/2) / (1-xso2*x/2); r1 =Po2*Pso2

4、/Pso3; r2 =Pso3/(Pso2*sqrt(Po2)*Kp); B =48148*exp(-7355.5/t); r3 =sqrt(B+(B-1)*(1-x)/x)+sqrt(K*(1-x)/x); r =keff*K*r1*(1-r2*r2)/(r3*r3); return(r);int main()double xso2=0.08,x,t,r0=0,y; int i; FILE *fp;fp=fopen("data.txt","w"); for(i=0;i<5;i+)fprintf(fp,"R=%en",r0);f

5、or(t=693.15;t<=873.15;)x=0.01;dox=x+h;if(x>=1) break;y=fabs(pow(10,5)*r(x,t,xso2)-pow(10,5)*r0);while(y>0.0001);fprintf(fp,"%5.2f %5.4fn",t,x);t=t+5;r0=r0+pow(10,-6);fclose(fp);return(0);運行上述程序后，獲得的數據見下表：溫度T/k反應速率0.0 1.0 ×10-62.0×10-6 3.0 ×10-64.0 ×10-6693.150.

6、9809 698.150.9797 703.150.9782 0.0117 708.150.9762 0.2719 713.150.9740 0.4750 718.150.9715 0.6205 0.1844 723.150.9686 0.7214 0.3863 0.0784 728.150.9655 0.7901 0.5373 0.2903 0.0638 733.150.9620 0.8365 0.6473 0.4536 0.2689 738.150.9582 0.8674 0.7265 0.5761 0.4285 743.150.9541 0.8873 0.7828 0.6669 0.55

7、00 748.150.9496 0.8995 0.8224 0.7336 0.6416 753.150.9447 0.8964 0.8237 0.7398 0.6529 758.150.9394 0.8927 0.8239 0.7444 0.6617 763.150.9338 0.8885 0.8231 0.7474 0.6685 768.150.9278 0.8837 0.8214 0.7491 0.6736 773.150.9213 0.8784 0.8188 0.7497 0.6772 778.150.9144 0.8725 0.8154 0.7492 0.6794 783.150.90

8、71 0.8662 0.8113 0.7477 0.6804 788.150.8994 0.8593 0.8065 0.7453 0.6803 793.150.8912 0.8519 0.8010 0.7421 0.6792 798.150.8826 0.8440 0.7948 0.7380 0.6772 803.150.8735 0.8357 0.7881 0.7332 0.6743 808.150.8640 0.8268 0.7807 0.7277 0.6706 813.150.8540 0.8175 0.7727 0.7215 0.6661 818.150.8436 0.8077 0.7

9、643 0.7146 0.6609 823.150.8327 0.7975 0.7552 0.7071 0.6550 828.150.8215 0.7868 0.7457 0.6990 0.6484 833.150.8097 0.7757 0.7357 0.6904 0.6412 838.150.7976 0.7642 0.7252 0.6812 0.6335 843.150.7851 0.7523 0.7142 0.6716 0.6252 848.150.7722 0.7400 0.7029 0.6614 0.6163 853.150.7590 0.7273 0.6912 0.6508 0.

10、6070 858.150.7454 0.7143 0.6791 0.6398 0.5972 863.150.7315 0.7010 0.6666 0.6285 0.5870 868.150.7174 0.6874 0.6539 0.6167 0.5764 873.150.7030 0.6736 0.6408 0.6047 0.5654 根據上述數據制圖： 上圖以溫度為橫坐標，轉化率為縱坐標。5個速率，0（對比，即反應平衡時）、1×10-6、2×10-6、3×10-6、4×10-6，作T-X的等速率線。因為是一維擬均相平推流，則；根據入口組成，設定入口溫度

11、，根據反應對入口溫度所求偏導數在這一段內對組成的積分為零，可以求得此段出口轉化率，即下一段的入口組成。催化劑用量最少時應使得；由于段內操作線的斜率為，因此根據入口溫度、入口組成和出口組成計算可得出口溫度。循環計算四段后可得最后的出口轉化率，若不能滿足要求（x98），則重新設定第一段入口溫度，再進行計算直至滿足條件。編寫c+程序，如下：#include<stdio.h>#include<math.h>#define R 1.987#define h 0.0001double r(double x, double t, double xso2)double r,r1,r2,

12、r3,keff,K,B,Kp,Pso2,Pso3,Po2; if(t>=693.15 && t<748.15) keff = 7.6915*pow(10,18)*exp(-76062 / (R*t); if(t>=748.15 && t<=873.15) keff = 1.5128*pow(10,7)*exp(-35992 / (R*t); K = 2.3*pow(10,-8)*exp(27200/(R*t); Kp=2.26203*pow(10,-5)*exp(11295.3 / t); Pso2=(xso2-xso2*x) / (1-

13、xso2*x/2); Pso3=xso2*x / (1-xso2*x/2);Po2=(0.17-xso2-xso2*x/2) / (1-xso2*x/2); r1 =Po2*Pso2/Pso3; r2 =Pso3/(Pso2*sqrt(Po2)*Kp); B =48148*exp(-7355.5/t); r3 =sqrt(B+(B-1)*(1-x)/x)+sqrt(K*(1-x)/x); r =keff*K*r1*(1-r2*r2)/(r3*r3); return(r);double dr(double x, double t, double xso2)double y; y=(r(x,t+h

14、,xso2)-r(x,t-h,xso2)/(2*h); return(y);double t(double t0, double x0, double x)double y,lamda,H=-23135,Cp=254.9,rou=0.500,c=1.282; lamda=-H*c/(rou*Cp); y=t0+lamda*(x-x0); return(y);double fun1(double x, double t, double xso2) double y; y=-dr(x,t,xso2)/(r(x,t,xso2)*r(x,t,xso2); return(y);double jifen(

15、double x0, double t0, double xso2) double sum=0.0,x1=x0,x2,t1,t2=693.15,xout; do t1=t(t0,x0,x1); x2=x1+h/10; t2=t(t0,x0,x2); if(t2>873.15) xout = x1; goto end; sum=sum+h*(fun1(x1,t1,xso2)+fun1(x2,t2,xso2)/20; x1=x2; while(sum<0);xout=x1-h/10;end: return(xout);double wjifen(double xin, double x

16、ou, double tin, double xso2)double y,x1=xin,x2,t1,t2,sum=0.0,wcat;do t1=t(tin,xin,x1); x2=x1+h; t2=t(tin,xin,x2); if(t2>=873.15) goto end ; sum=sum+(1/r(x1,t1,xso2)+1/r(x2,t2,xso2)*h/1000; x1=x2; while(x2<=xou);end: wcat=sum*131*1000/3600; return(wcat);int main()double xso2=0.08,xout,tout,x0=0

17、.0001,t0,t00=719,t1,wcat,wsum;int i,j;loop1:wsum=0.0;x0=0.0001;t0=t00;FILE*fp;fp=fopen("data2.txt","w");fprintf(fp,"1 tin=%f xin=%en",t00,x0);for(i=0;i<=3;i+)xout=jifen(x0,t0,xso2);tout=t(t0,x0,xout);fprintf(fp,"%d tout=%f ",i+1,tout);fprintf(fp,"xout=

18、%f n",xout);wcat=wjifen(x0,xout,t0,xso2);fprintf(fp,"Wcat=%f ",wcat);wsum=wsum+wcat;t1= 693.15;dot1=t1+0.01;while(fabs(pow(10,5)*r(xout,t1,xso2)-pow(10,5)*r(xout,tout,xso2)>h);x0=xout;t0=t1;fprintf(fp,"%d tin= %f xin= %fn",i+2,t0,xout);fprintf(fp,"Wsum=%f n",wsu

19、m);t00=t00-0.1;fprintf(fp,"nn");if(x0<=0.98) goto loop1;return(0);所得數據如下：第一節反應器： tin=717.800000 xin=1.000000e-004 tout=873.148833 xout=0.667660 Wcat=8349.667114（kg） 第二節反應器： tin= 723.600000 xin= 0.667660 tout=778.519894 xout=0.903660 Wcat=11605.588273（kg） 第三節反應器： tin= 715.920000 xin= 0.9

20、03660 tout=729.505692 xout=0.962040 Wcat=21557.219801（kg） 第四節反應器： tin= 693.190000 xin= 0.962040 tout=697.388114 xout=0.980080 Wcat=50379.460102（kg） 總催化劑：Wsum=91891.935289（kg） 作折線圖如下：討論：（1） 將第二個程序改為99%，因為在方程給出的最低溫度693.15K時轉化率為98.38%，無法達到99%，所以得不出99%轉化率時的催化劑用量。（2） 程序如下（舉例，YSO2=7%）：#include"math.h

21、"#include"stdio.h"#define R 1.987#define h 0.0001double r(double x,double t,double xso2)double r,r1,r2,r3,keff,K,B,Kp,Pso2,Pso3,Po2; if(t>=693.15 && t<748.15) keff = 7.6915*pow(10,18)*exp(-76062 / (R*t); if(t>=748.15 && t<=873.15) keff = 1.5128*pow(10,7)*ex

22、p(-35992 / (R*t); K = 2.3*pow(10,-8)*exp(27200/(R*t); Kp=2.26203*pow(10,-5)*exp(11295.3 / t); Pso2=(xso2-xso2*x) / (1-xso2*x/2); Pso3=xso2*x / (1-xso2*x/2); Po2=(0.21-xso2-xso2*x/2) / (1-xso2*x/2); r1 =Po2*Pso2/Pso3; r2 =Pso3/(Pso2*sqrt(Po2)*Kp); B =48148*exp(-7355.5/t); r3 =sqrt(B+(B-1)*(1-x)/x)+sq

23、rt(K*(1-x)/x); r =keff*K*r1*(1-r2*r2)/(r3*r3); return(r);double dr(double x,double t,double xso2)double y; y=(r(x,t+h,xso2)-r(x,t-h,xso2)/(2*h); return(y);double t(double t0,double x0,double x)double y,lamda,H=-23135,Cp=254.9,rou=0.500,c=1.282; lamda=-H*c/(rou*Cp); y=t0+lamda*(x-x0); return(y);doubl

24、e fun1(double x,double t,double xso2) double y; y=-dr(x,t,xso2)/(r(x,t,xso2)*r(x,t,xso2); return(y);double jifen(double x0,double t0,double xso2) double sum=0.0,x1=x0,x2,t1,t2=693.15,xout; do t1=t(t0,x0,x1); x2=x1+h/10; t2=t(t0,x0,x2); if(t2>873.15) xout = x1; goto end; sum=sum+h*(fun1(x1,t1,xso2

25、)+fun1(x2,t2,xso2)/20; x1=x2; while(sum<0);xout=x1-h/10;end: return(xout);double wjifen(double xin,double xou,double tin,double xso2)double y,x1=xin,x2,t1,t2,sum=0.0,wcat;do t1=t(tin,xin,x1); x2=x1+h; t2=t(tin,xin,x2); if(t2>=873.15) goto end ; sum=sum+(1/r(x1,t1,xso2)+1/r(x2,t2,xso2)*h/1000;

26、x1=x2; while(x2<=xou);end: wcat=sum*131*1000/3600; return(wcat);int main() double xso2=0.07,xout,tout,x0=0.0001,t0,t00=717,t1,wsum=0.0,wcat; double xou5,tou5,tin5,xin5; int i,j; FILE*fp; fp=fopen("data4.txt","w");loop1: x0=0.0001; t0=t00; i=0; xini=x0; tini=t00;loop2: xout=jifen(x0,t0,xso2); tout=t(t0,x0,xout); xoui=xout; toui=tout;t1= 693.15;do t1=t1+0.01; while(fabs(pow(10,5)*r(xout,t1,xso2)-pow(10,5)*r(xout,tout,xso2)>h);x0=xout; t0=t1;i+; xini=xout; tini=t0; if(i<=3) goto loop2;t00=t00-0.1;if(x0<=0.98) goto loop1; fprintf(fp,"