SO2反应器最优化(北化化反大作业)

ww--二氧化硫转化器最优化声明Fortran90matlabCFortranFortranmatlabCSO1/2O=SO2 2 3Cp0.2549[kcal/k·-ΔH=23135[kcal/kmol]bρm]bS2

9.0。131[kmolS28。 动力 －R k P

KPSO2

PSO3

1

/其

SO2

eff

BB1PSO2

P SO3

KPSO2

P 2SO376062 keff

7.6915exp

RT

420475oC35992 keff

1.5128107exp

RT

475600oCB48148exp7355.5 T 27200K2.3108exp

RT P KP

P12SO P3

SO2

O211295.3K 2.26203105exp P T

R1.987T~X图上做出平衡线4条等速率线平衡线：A1/2BCr02=01A整理： PC =KpPAPB1/2其中： eXA0~1范围内变化XAe唯一对的TXA从ee1.0开始间隔-0.1分别算出相T即线：e等速率线

图1 温度—转化率示意图rX

能对TTA AX绘制等速率线r再A ．．．．A．A。然A也可． ．。但是此处绘制速率多的是为了得到的趋势而的降低并不影响曲的走势。图2 ~0℃ 线要求 2：以一维拟均相平推流模型为基础，在催化剂用量最少的前提下，总的及各段的催化剂用量；进出口温度、转化率；并在 T 图上标出折线。分析：要使总得催化剂用量最少，则每一段催化剂用量都要达到最小。当确定每一段的Xin

和Xout

之后，根据0

可确定最佳进口温度 T,

，使得 ,最小，即满足 ：

0（式 1）。但在实际计算中，第一段的 in是知道的，只需假定一个Tin值，即可根据上式得到 t，进而由操作线算出Tout。而段与段之间有

，为使任意两段之间催化剂之和最小，进一, ,步可以得出段与段之间需满足： , , 由, , , ,上式可确定下一段入口温度Tin,i+1，如此循环四次即可得到各段反应器进出口温度、转化率，以及催化剂用量。．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．A．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．’．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．以下几计算中的问题，应该得到重视：XinTin计算Xout有两个限制条件，一是能越Xout和操作线确定的温度能越600℃②根据 , , ,

,, ,

求解 Tin,i+1时，使用割线迭代初值要求很严格，如果初值选择合适，甚至无解。虽然有其他高阶算可以避免处选择的是一个“笨办，可以叫做试探：T0℃开始，微小增量进行迭代计算，直到满足 ri+1=ri为止。此虽然看起来“笨”，精度低，但却在一定度上加快序运行速度，巧妙得解决问题。~0T1Xt8Tn1最终转化率与入口温度关系T℃ Xin out4410.980094420.980244430.980294440.98027~5小T5Xout8nXAXAWcat

4750047000XA0.98015

Wcat0.98010

450000.980050.98000

44500440 441 442 443 444 T /℃in图3~从上图看出98%5%44541.6Kg2X T/℃

T /℃

/Kg

/Kgin

out

out

cat

cat10.0000010.00000444.550.66917600.004146.7720.66917450.360.90470505.085767.1530.90470442.740.96243456.1410670.2440.96243420.380.98007424.4823957.471.00.80.6XA0.40.20.0420 440 460 480 500 520 540 560 580 T/℃图4 9999%99Te0~099%。2 2 O＋O＝1O进口浓－9之间装改02 2 。3浓与总关系0 Tn℃ cat0.07465.4012224.200.08443.7516107.800.09423.7523355.99＋SO＝21以及达98前提随着入口SO2浓增加总也增加。源代码：计算过程中所有子程序ModuleCR_1ImplicitNoneReal(kind=8)::PA,PB,PC,Kp,T,keff,B,K,ksi,rA_C,&rA,DT,DX,FT,Xmax,S,Xin,Xout,Tin,Tout,&XA,W,X,TOL,XAe,N,E,FX,det,X0,FX0,&Sum,W_sum,Tin_1,lmdReal,parameter::P=1.01325,yA0=0.08,yB0=0.09,R=1.987Contains子程序：计算各组分分压SubroutineCal_P()PA=yA0*P*(1-XA)/(1-0.5*yA0*XA)PB=P*(yB0-0.5*yA0*XA)/(1-0.5*yA0*XA)PC=yA0*P*XA/(1-0.5*yA0*XA)EndSubroutineCal_P子程序：计算反应速率SubroutineCal_rA()IF(T<748.15)Thenkeff=7.6915E18*EXP(-76062/R/T)Elsekeff=1.5128E7*EXP(-35992/R/T)EndIFB=48148*EXP(-7355.5/T)K=2.3E-8*EXP(27200/R/T)Kp=2.26203E-5*EXP(11295.3/T)ksi=PC/Kp/PA/SQRT(PB)rA_C=keff*PB*K*PA/PC*(1-ksi**2)/(SQRT(B+(B-1)*PA/PC)+SQRT(K*PA/PC))**2EndSubroutineCal_rASubroutineSecant()IF(N<=1)Thendet=-0.1Elsedet=-FX*(X-X0)/(FX-FX0)EndIFX0=XFX0=FXE=abs(det)/(abs(X)+1.)X=X+detEndSubroutineSecantSubroutineCal_XAe()X=1N=0130 N=N+1FX=2.26203E-5*exp(11295.3/(lmd*(X-Xin)+Tin))*(1-X)*sqrt(P*(yB0-yA0*X/2)/&(1-yA0*X/2))-XCallSecant()IF(E>=1.E-6)GOTO130XAe=XEndSubroutineCal_XAe一段由T，XXin SubroutineCal_Xout()S=0DT=0.000001DX=0.00001T=TinXA=XinCallCal_XAe()

outXmax=(873.15-T)/lmd+XinIF(Xmax>1.0)Xmax=1.0100 CallCal_P()CallCal_rA()rA=rA_CT=T+DTCallCal_rA()FT=(rA_C-rA)/DT/rA**2S=S+FT*DXIF(XA>=XAe)ReturnIF(XA>=Xmax)ReturnIF(S>0)ThenXA=XA+DXT=lmd*DX+TGOTO100EndIFEndSubroutineCal_Xout一段催化剂装量SubroutineCal_W()Sum=0T=TinXA=XinDX=0.00001110 CallCal_P()CallCal_rA()rA=rA_CFT=1/rASum=Sum+FT*DXIF(XA>=Xout)ReturnXA=XA+DXT=T+lmd*DXGOTO110EndSubroutineCal_W由上一段T，X一段ToutSubroutineCal_T()XA=XoutT=ToutCallCal_P()CallCal_rA()rA=rA_CT=420+273.15120 CallCal_rA()

out inIF(abs(rA-rA_C)<1.E-11)ReturnEndIFT=T+0.0001GOTO120EndSubroutineCal_TEndModuleCR_1主调用程序：主程序：计算平衡线ModuleCR_0ImplicitNoneReal(kind=8) ::PA,PB,PC ,Kp,TReal::XAReal,parameter::P=1.01325,yA0=0.08,yB0=0.09ContainsSubroutineCal_1()PA=yA0*P*(1-XA)/(1-0.5*yA0*XA)PB=(yB0*P-0.5*yA0*P*XA)/(1-0.5*yA0*XA)PC=yA0*P*XA/(1-0.5*yA0*XA)Kp=PC/PA/SQRT(PB)T=11295.3/LOG(Kp*100000/2.26203)EndSubroutineCal_1EndModuleCR_1ProgramTe-XAeuseCR_0ImplicitNoneOpen(13,file='data1.csv')XA=0.05100 DoWhile(XA<=1.0)CallCal_1()IF(T>693.15.AND.T<873.15)ThenWrite(13,*)XA,',',T-273.15IF(XA<0.95)ThenXA=XA+0.05ElseXA=XA+0.001EndIFGOTO100EndDoEndProgramTe-XAe主程序：等速率曲线ProgramrA_XA_TuseCR_1ImplicitNoneInteger :: Real::rAOpen(12,file='data2.csv')Doi=1,5rA=i*1.e-6Write(12,*)"rA=",rAT=693.15200XA=0.001100CallCal_P()CallCal_rA()IF(ABS(rA_C-rA)<1.E-10)ThenWrite(12,*)"XA",',',XA,',',"T",',',T-273.15EndIFXA=XA+0.00001IF(XA<1.0)GOTO100T=T+5IF(T<=883.15)GOTO200EndDoEndProgramrA_XA_T主程序：输出所有转化率>98%的结果ProgrammainuseCR_1ImplicitNoneInteger :: Open(15,file='Tin-Xout.csv')lmd=23135*yA0/0.2549/(64.065*yA0+31.999*yB0+(1-yA0-yB0)*28.013)Tin_1=420+273.15300Tin=Tin_1Xin=1.E-10W_sum=0i=0200i=i+1CallCal_Xout()Xout=XATout=lmd*(Xout-Xin)+TinIF(Tout>873.15)Tout=873.14999CallCal_W()W=Sum*131/3600W_sum=W_sum+WCallCal_T()IF(abs(T-Tout)<0.5)ThenIF(i<4)GOTO400EndIFTin=TXin=XoutIF(i<4)GOTO200IF(Xout>0.98)ThenWrite(15,*)Tin_1-273.15,',',Xout,',',W_sumE