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

《化学反应器理论》大作业化研1512班xxx2015200xxx1二氧化硫转化器最优化声明：本人使用Fortran90语言独立完成该作业。在前几届学生遗留下来的作业模板中，多数使用matlab软件或C语言编译，本人至今未见使用Fortran语言编译的作业。私以为Fortran语言相比于matlab语言或C语言，有着得天独厚的优势，更具有学习价值。故本作业的独立性毋庸置疑！题目：SO2＋1/2O2=SO3，四段绝热反应器，级间间接换热。1.基础数据：混合物恒压热容Cp＝0.2549[kcal/kg·K]-ΔH=23135[kcal/kmol]床层空隙率ρb＝554[kg/m3]进口SO2浓度8.0%，O2浓度9.0%，其余为氮气。处理量131[kmolSO2/hr]，要求最终转化率98％。2.动力学方程：其中：sec./112232323222gcatmolPPKPPBBPPKPkRSOSOSOSOSOSOOeffSO－987.13.11295exp1026203.227200exp103.25.7355exp4814860047535992exp105128.147542076062exp106915.75218718223RTKPPKPRTKTBCRTkCRTkPOSOPSOoeffoeff《化学反应器理论》大作业化研1512班xxx2015200xxx2要求1：在T~X图上，做出平衡线，至少4条等速率线（1）平衡线：反应简写为：A+1/2B=C当反应平衡时，（-rA）=0，得到（1-ξ2）=0，即ξ=1整理得：PC=KpPAPB1/2其中：00.08Ay00.09By0.5A当XA从0~1范围内变化时，可确定与XAe唯一对应的平衡温度Te，编写程序使XA从1.0开始，间隔为-0.1，分别算出相应的Te，即得平衡线：图1平衡温度—平衡转化率示意图00(1)1AAAAAAyPXPyX0001(1)21BAABAAAyPyPXPyX001AACAAAyPXPyX《化学反应器理论》大作业化研1512班xxx2015200xxx3（2）等速率线由于在同一反应速率（-rA）下，某一转化率XA可能对应两个温度T，而某一温度T下，则对应唯一转化率XA。所以．．绘制等速率线时，．．．．．．．．我．采用的方法为：．．．．．．．先．确定反．．．应速率（．．．．-．r．A．）．，再．．确定．．温度．．T．，然后对．．．．X．A．从．0．到．1．进行迭代，．．．．．找到．．那个唯一的．．．．．X．A．。当然，在求解XA时也可以使用割线法等其他迭代算法，那样精度更高。但是此处绘制等速率线更多的是为了得到XA~T的趋势，而精度的降低并不影响曲线的走势。图2420~600℃下平衡线与等反应速率线《化学反应器理论》大作业化研1512班xxx2015200xxx4要求2：以一维拟均相平推流模型为基础，在催化剂用量最少的前提下，总的及各段的催化剂用量；进出口温度、转化率；并在T-X图上标出折线。分析：要使总得催化剂用量最少，则每一段催化剂用量都要达到最小。当确定每一段的Xin和Xout之后，根据0(,)outinXcatAXAAAWdXFrXT可确定最佳进口温度Tin，使得Wcat最小，即满足：(,)10outinXAAAXAinrXTdXrT（式1）。但在实际计算中，第一段的Xin是知道的，只需假定一个Tin值，即可根据上式得到Xout，进而由操作线算出out。而段与段之间有,,1outiiniXX，为使任意两段之间催化剂之和最小，进一步可以得出段与段之间需满足：,,,1,1(,)(,)outioutiiniinirXTrXT由上式可确定下一段入口温度Tin,i+1，如此循环四次即可得到各段反应器进出口温度、转化率，以及催化剂用量。需要特别注意的是：．．．．．．．．．在．计算积分式时，．．．．．．．采用．．的是数值微积分方法。．．．．．．．．．．在对．．X．A．取微小步．．．．长．Δ．X．时．，．T．也会沿着操作线移动．．．．．．．．．，．得到．．新的．．T．’．，．下一步．．．的数值微积分是基于新的．．．．．．．．．．．T．’．而计算的．．．．，．如此．．循环往复．．．．。．在．编程的过程中．．．．．．，．正．是由于本人忽略了这一点．．．．．．．．．．．，．误．认为．．T．值是不变的．．．．．，．从而导致了严重的计算错误．．．．．．．．．．．．，．始终．．得不到正确的解．．．．．．．，．浪费了．．．大量时间查．．．．．找算法上的错误．．．．．．．，．而．忽略了对问题本质的研究．．．．．．．．．．．。．这也．．让我想起了导师王健红的一句话．．．．．．．．．．．．．．：．“．建模．．的关键不是优化算法，而是准确认识模型的原型！．．．．．．．．．．．．．．．．．．．．．．”．以下几点计算过程中的问题，也应该得到重视：①由Xin和Tin计算Xout时,Xout有两个限制条件，一是不能越过平衡线，二是Xout和操作线确定的温度不能越过600℃②根据,,,1,1(,)(,)outioutiiniinirXTrXT求解Tin,i+1时，使用割线法对迭代初值要求很严格，如果初值选择不合适，会导致得不到想要的解，甚至无解。虽然有其他高阶算法可以避免这类问题，但由于反应速率方程难以求得解析导，故不适用。此处本人选择的是一个“笨办法”，也可以叫做试探法：T从420℃开始，取微小增量进行迭代计算，直到满足ri+1=ri为止。此法虽然看起来“笨”，精度略低，但却在一定程度上加快了程序运行速度，巧妙得解决了问题。《化学反应器理论》大作业化研1512班xxx2015200xxx5解题思路：第一步，在420~600℃大范围内，以ΔT=1逐步计算，输出所有Xout0.98的Tin值。程序运行结果如下表所示。表1最终转化率与入口温度关系Tin/℃Xout4410.980094420.980244430.980294440.98027第二步，在440~445℃小范围内，以ΔT=0.05逐步计算，输出所有Xout0.98的Tin值。程序运行结果如下图所示。4404414424434444450.980000.980050.980100.980150.980200.980250.980300.980350.98040XAWcatTin/℃XA44500450004550046000465004700047500Wcat图3440~445℃范围内最终转化率及催化剂总量与入口温度的关系从上图可以看出，在满足最终转化率达到98%的前提下，入口温度越高，催化剂总用量越少，因此选择最佳入口温度为444.55℃，此时出口转化率为98.01%，催化剂总用量为44541.6Kg。各段反应器具体操作指标如下表所示。《化学反应器理论》大作业化研1512班xxx2015200xxx6表2最终转化率达标、催化剂总量最少条件下各段反应器指标XinTin/℃XoutTout/℃Wcat/KgWcat/Kg10.00000444.550.66917600.004146.7744541.6220.66917450.360.90470505.085767.1530.90470442.740.96243456.1410670.2440.96243420.380.98007424.4823957.474204404604805005205405605806000.00.20.40.60.81.0XAT/℃图4T~XA图中操作线《化学反应器理论》大作业化研1512班xxx2015200xxx7讨论1：要求的最终转化率从98％变化到99％对催化剂用量的影响从程序运行的角度看，始终没有得到转化率超过99%的解。而从问题的本质来看，当转化率为99%时，平衡温度Te为400℃。而催化剂的适用温度为420~600℃，故该题最终转化率无法达到99%。讨论2：yO2＋ySO2＝21％，SO2进口浓度在7－9％之间变化，对催化剂装量的影响改变yA0的初值，运行程序，输出所有满足转化率98%的结果，并比较其催化剂总量，找出催化剂用量最少的结果，并记录对应的入口温度值，得到下表。表3入口SO2浓度与催化剂总量的关系yA0Tin/℃Wcat/Kg0.07465.4012224.200.08443.7516107.800.09423.7523355.99从上表可以看出，在满足yO2＋ySO2＝21％以及最终转化率达到98%的前提下，随着入口SO2浓度的增加，催化剂总量也增加。《化学反应器理论》大作业化研1512班xxx2015200xxx8源代码：计算过程中所有子程序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(T748.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_rA子程序：割线法求根SubroutineSecant()IF(N=1)Thendet=-0.1Elsedet=-FX*(X-X0)/(FX-FX0)EndIFX0=XFX0=FXE=abs(det)/(abs(X)+1.)X=X+detEndSubroutineSecant《化学反应器理论》大作业化研1512班xxx2015200xxx9子程序：计算操作线与平衡线交点的平衡转化率SubroutineCal_XAe()X=1N=0130N=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子程序：每一段由Tin，Xin计算XoutSubroutineCal_Xout()S=0DT=0.000001DX=0.00001T=TinXA=XinCallCal_XAe()Xmax=(873.15-T)/lmd+XinIF(Xmax1.0)Xmax=1.0100CallCal_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(S0)ThenXA=XA+DXT=lmd*DX+TGOTO100EndI