反应器大作业(SO2优化)
化学反应器理论大作业——二氧化硫转化器最优化
学院：化学工程学院
班级：化研1108
学号：2011200126
姓名：马树刚
二氧化硫转化器的最优化
题目背景：SO 2＋1/2O 2=SO 3，四段绝热反应器，级间间接换热。
1. 基础数据：
混合物恒压热容 C p ＝0.2549[kcal/kg·K]
－ΔH =23135[kcal/kmol]
床层空隙率 ρb ＝554[kg/m 3]
进口SO 2浓度8.0%，O 2浓度9.0%，其余为氮气。 处理量131[kmolSO 2/hr]，要求最终转化率98％。
2. 动力学方程：
式中：
3. 基本要求：
(1)在T －X 图上，做出平衡线，至少4条等速率线；
(2)以一维拟均相平推流模型为基础，在催化剂用量最少的前提下，总的及各段的催化剂装量；进出口温度、转化率；并在T-X 图上标出折线； (3)程序用C ，Fortran ，BASIC 语言之一编制；
()
()(
)
[]
sec ./112
2
3
232322
2gcat mol P P K P P B B P P K P k R SO SO SO SO SO SO O eff SO +-+-=ξ－()
()
987
.13.11295
exp 1026203.227200exp 103.25.7355exp 4814860047535992exp 105128.147542076062exp 106915.7521
8718223
=?
?
? ???=??
? ??=?
?
?
???=?
??
??-=-?
??
??-?=-?
??
??-?=--R T K P P K P RT K T B C RT k C RT k P O SO P SO
o
eff o
eff ξ
4.讨论：
(1)要求的最终转化率从97％变化到99％对催化剂用量的影响；
(2)如果有关系：YO2＋YSO2＝21％，SO2进口浓度在7－9％之间变化，对催化
剂装量的影响。
5.选做：对优化结果进行校核。
算法简介及计算程序和结果
1.解题思路
已知第一段入口和最后一段出口的转化率；第一段入口反应物浓度，各物性参数；段与段间采用间接冷却。可改变的参数为各段的入口温度；段与段之间的转化率。
由上述反应流程图分析可知：根据已知的入口组成，设定入口温度，根据反应速率对入口温度所求偏导数在这一段内对组成的积分为零可以求得此段出口转化率和出口温度，即得到下一段的入口转化率，又根据前一段的出口速率等于后一段的入口速率，可以求得下一段的入口温度；这样又可以计算下一段的出口情况。这样反复计算直至求出最后的出口情况。判断此出口情况是否满足题给条件（转化率达到98%），若不满足条件，则应重新假设再进行计算。这样反复寻优直至满足条件。计算公式如下：
第一段：
0111
2
111=???
dx T r r in
X X out
in
第一、二段之间：),(),(222111in in out out T X r T X r = 第二段：
0122
2
222=???
dx T r r in
X X out
in
第二、三段之间：),(),(333222in in out out T X r T X r = 第三段：
0133
2333=???
dx T r r in
X X out
in
第三、四段之间：),(),(444333in in out out T X r T X r = 第四段：
0144
2
444=???
dx T r r in
X X out
in
)(O O X X T T -=-λ 其中：λ为绝热温升
2. 计算程序和结果
2.1平衡线及等速率线的计算程序：
利用反应速度表达式，代入不同的速度值来确定不同温度下的转化率，从而求取等速率线和平衡线（速度等于0）。程序如下：
#include "math.h" #include "stdio.h" #include "iostream.h" #include "stdlib.h" const double R=1.987;
void strcopy(char *str1,char *str2) {
for(int i=0; i<6; i++) { str1[i]=str2[i]; } }
double rate( double x, double t) {
double y,y1,y2,y3,keff,k,B; 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));
y1 = ((0.09-0.04*x) / (1-0.04*x)) * ((0.08-0.08*x) / (1-0.04*x)) / (0.08*x / (1-0.04*x));
y2 = (0.08*x / (1-0.04*x)) / (((0.08-0.08*x) / (1-0.04*x)) * sqrt((0.09-0.04*x) / (1-0.04*x)) * 2.26203*pow(10,-5)*exp(11295.3 / t));
B = 48148*exp(-7355.5 / t);
y3 = sqrt(B+(B-1)*(1-x) / x) + sqrt(k*(1-x) / x);
y = keff*k*y1*(1-y2*y2) / (y3*y3);
return y;
}
void main()
{
double x,t,d;
double r0=0,step=0.0001;
FILE *fp;
char filename[6];
for(int i=0; i<=5; i++)
{
t=693.15;
do
{
x=0.01;
do
{
x+=step;
if(x>=1)
goto next;
d=fabs(pow(10,5)*rate(x,t) - pow(10,5)*r0);
}
while(d>0.001);
switch(i)
{
case 0:
strcopy(filename,"data0");
break;
case 1:
strcopy(filename,"data1");
break;
case 2:
strcopy(filename,"data2");
break;
case 3:
strcopy(filename,"data3");
break;
case 4:
strcopy(filename,"data4");
break;
default:
strcopy(filename,"data5");
}
fp=fopen(filename,"a");
if ((x-0.0001)<1)
{
printf("Xso2= %f T= %f\n",x-0.0001,t);
fprintf(fp,"%f %f\n",x-0.0001,t);
fclose(fp);
}
next: t+=5;
}
while(t<=873.15);
r0 +=1.0*pow(10,-6);
}
}
2.2 催化剂用量最小时的温度、转化率及催化剂装量计算程序
在催化剂用量最少（又能满足要求的最高反应速率）的前提下，计算各段进出
口温度、转化率以及总的和各段的催化剂装量。程序如下：
#include"math.h"
#include"stdio.h"
#define R 1.987
#define h 0.0001
double myabs(double x)
{
return(x>0?x:(-x));
}
double r(double x,double t)
{
double y,y1,y2,y3,keff,k,b;
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));
y1 = ((0.09-0.04*x) / (1-0.04*x)) * ((0.08-0.08*x) / (1-0.04*x)) / (0.08*x / (1-0.04*x));
y2 = (0.08*x / (1-0.04*x)) / (((0.08-0.08*x) / (1-0.04*x)) * sqrt((0.09-0.04*x) / (1-0.04*x)) * 2.26203*pow(10,-5)*exp(11295.3 / t));
b = 48148*exp(-7355.5 / t);
y3 = sqrt(b+(b-1)*(1-x) / x) + sqrt(k*(1-x) / x);
y = keff*k*y1*(1-y2*y2) / (y3*y3);
return(y);
}
double dr(double x,double t)
{
double y;
y= (r(x,t+h) - r(x,t-h)) / (2*h);
return(y);
}
double t(double t0,double x0,double x)
{
double y,lamda,H,cp,rou,c;
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 y;
y= -dr(x,t) / (r(x,t)*r(x,t));
return(y);
}
double jifen(double x0,double t0)
{
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) + fun1(x2,t2)) / 20;
x1=x2;
}
while(sum<0);
xout = x1-h/10;
end: return(xout);
}
double wjifen(double xin,double xou,double tin)
{
double x1=xin,x2,t1,t2,sum=0.0,wcat;
do
{
t1=t(tin,xin,x1);
x2=x1+h;
t2=t(tin,xin,x2);
sum = sum + (1/r(x1,t1) + 1/r(x2,t2))*h / 2000;
x1=x2;
}
while(x2<=xou);
wcat = sum*131*1000/3600;
return(wcat);
}
void main()
{
double xout,tout,x0=0.0001,t0,t00=720,t1,wsum=0.0,wcat;
double xou[5],tin[5],xin[5],tou[5];
int i,j;
do
{
x0 = 0.0001;
t0 = t00;
printf("1 tin= %f xin= %e\n",t00,x0);
i=0;
xin[i] =x0;
tin[i] =t00;
do
{
xout = jifen(x0,t0);
tout = t(t0,x0,xout);
xou[i] = xout;
tou[i]=tout;
printf("%d tout= %f ",i+1,tout);
printf("xout= %e \n",xout);
printf("*****************************************\n");
t1= 693.15;
do
{
t1 = t1 + 0.01;
}
while(myabs(pow(10,5)*r(xout,t1)-pow(10,5)*r(xout,tout))>h);
x0 = xout;
t0 = t1;
printf("%d tin= %f xin= %e\n",i+2,t0,xout);
i++;
xin[i] = xout;
tin[i] = t0;
/* printf("%d tin[]= %f xin[]= %e\n",i+1,tin[i],xin[i]); */
}
while(i<=3);
t00 = t00-0.1;
printf("\n\n");
}
while(x0<=0.98);
FILE *fp;
fp=fopen("data.txt","w");
for(j=0;j<=3;j++)
{
printf("%d %f %e %f %e\n",j+1,tin[j],xin[j],tou[j],xou[j]);
fprintf(fp,"%d tin=%f xin=%e tout=%f xou=%e\n",j+1,tin[j],xin[j],tou[j],xou[j]);
wcat = wjifen(xin[j],xou[j],tin[j]);
printf("%d wcat=%f\n",j+1,wcat);
fprintf(fp,"%d wcat=%f\n",j+1,wcat);
wsum = wsum + wcat;
}
printf("wsum= %f\n",wsum);
fclose(fp);
}
2.3 运行结果
(1)等速线和平衡线程序运算结果：（其中R为反应速率；T（K）为反应温度;
表内主体为各反应速率下对应系列温度的SO2转化率；）
表2.1
(2) 关于催化剂用量程序运算结果
以下为运行结果：（经运算，满足给定条件的第一段最高入口温度为717.80K，此时可得到最少的催化剂用量）
表2.2
三、平衡线、等速线及操作线图
根据以上数据利用Matlab进行编辑得图，程序如下：%平衡线和等速线
load data0;
x=data0(:,1);
t=data0(:,2);
plot(t,x,'r');
load data1;
x=data1(:,1);
t=data1(:,2);
hold on
plot(t,x,'g');
load data2;
x=data2(:,1);
t=data2(:,2);
hold on
plot(t,x);
load data3;
x=data3(:,1);
t=data3(:,2);
hold on
plot(t,x,'y');
load data4;
x=data4(:,1);
t=data4(:,2);
plot(t,x);
hold on
load data5;
x=data5(:,1);
t=data5(:,2);
hold on
plot(t,x);
%操作线
y=[10^(-4) 717.800000;
6.6766e-001 873.148833 ;
6.676600e-001 723.600000;
9.036600e-001 778.519894;
9.036600e-001 715.920000;
9.620400e-001 729.505692;
9.620400e-001 693.190000;
9.800800e-001 697.388114];
hold on
plot(y(:,2),y(:,1),'-.r');
得图如下：
图3.1平衡线、等速线及操作线图
四、讨论
1. 由计算结果可知转化率由97%变化到98%时,随着转化率的升高（98%变化到
99%的结果程序中无法输出，但是由分析可知其结论是一样），催化剂的用量显著升高，上升的幅度很大。即当转化率很高时，再提高很小的转化率时，结果所需的催化剂的量大大增加。这可能由于达到较高转化率时进一步提高转化率，反应物浓度很低，反应物被吸附到催化剂表面参与反应的几率降低，因而要继续提高转化率需保证更大催化接触表面，因而催化剂用量大幅增加。由如下计算结果即可知。 X=97.0%时，Weightcat=23.491T
X=97.5%时，Weightcat=37.890T X=98.0%时，
Weightcat=45.946T
图4.1 SO2转化率与催化剂用量关系图
2、YO2＋YSO2＝21％，SO2进口浓度在7－9％之间变化时：
YSO2=0.070时，Wcat=26.087T
YSO2=0.075时，Wcat=29.196T
YSO2=0.080时，Wcat=33.108T
YSO2=0.085时，Wcat=38.125T
Wcat=45.467T
YSO2=0.090时，
可以看出，随着SO2入口含量的增加，催化剂的用量也随之增加，而且其变化的相对量也在增加，其原因可能是由于在SO2和O2总含量不变的情况下，因为O2相对SO2大大过量，增加SO2含量相当于提高了反应物浓度。在转化率不变条件下，需要催化剂提供更大表面积作为它们的反应场所，故而需要更多的催化剂。