SO2反應器最優(yōu)化(北化化反大作業(yè))

ww--二氧化硫轉(zhuǎn)化器最優(yōu)化聲明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范圍內(nèi)變化XAe唯一對的TXA從ee1.0開始間隔-0.1分別算出相T即線：e等速率線

圖1 溫度—轉(zhuǎn)化率示意圖rX

能對TTA AX繪制等速率線r再A ．．．．A．A。然A也可． ．。但是此處繪制速率多的是為了得到的趨勢而的降低并不影響曲的走勢。圖2 ~0℃ 線要求 2：以一維擬均相平推流模型為基礎(chǔ)，在催化劑用量最少的前提下，總的及各段的催化劑用量；進出口溫度、轉(zhuǎn)化率；并在 T 圖上標出折線。分析：要使總得催化劑用量最少，則每一段催化劑用量都要達到最小。當確定每一段的Xin

和Xout

之后，根據(jù)0

可確定最佳進口溫度 T,

，使得 ,最小，即滿足 ：

0（式 1）。但在實際計算中，第一段的 in是知道的，只需假定一個Tin值，即可根據(jù)上式得到 t，進而由操作線算出Tout。而段與段之間有

，為使任意兩段之間催化劑之和最小，進一, ,步可以得出段與段之間需滿足： , , 由, , , ,上式可確定下一段入口溫度Tin,i+1，如此循環(huán)四次即可得到各段反應器進出口溫度、轉(zhuǎn)化率，以及催化劑用量。．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．A．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．’．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．以下幾計算中的問題，應該得到重視：XinTin計算Xout有兩個限制條件，一是能越Xout和操作線確定的溫度能越600℃②根據(jù) , , ,

,, ,

求解 Tin,i+1時，使用割線迭代初值要求很嚴格，如果初值選擇合適，甚至無解。雖然有其他高階算可以避免處選擇的是一個“笨辦，可以叫做試探：T0℃開始，微小增量進行迭代計算，直到滿足 ri+1=ri為止。此雖然看起來“笨”，精度低，但卻在一定度上加快序運行速度，巧妙得解決問題。~0T1Xt8Tn1最終轉(zhuǎn)化率與入口溫度關(guān)系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濃與總關(guān)系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主調(diào)用程序：主程序：計算平衡線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主程序：輸出所有轉(zhuǎn)化率>98%的結(jié)果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