1、化學反應器理論大作業二氧化硫轉換器最優化班級 :化研 1612學號 :2016200162指導教師:文利雄姓名 : 閆曉宇化學反應器理論大作業二氧化硫轉換器最優化反應方程式：SO2 1/2O2=SO3 （放熱反應）四段絕熱反應器，級間間接換熱，常壓下反應。1. 基礎數據? 混合物恒壓熱容Cp 0.2549 kcal/kgK? H =23135 kcal/kmol? 催化劑堆密度 b 554 kg/m 3? 進口SO2濃度8.0 mol%， O2濃度9.0mol%，其余為氮氣? 處理量 131 kmol SO 2/hr ，要求最終轉化率98KPSO2 PSO3 12RSO2keff PO2mol

2、/ gcat. sec1.1. 動力學方程2B B 1 PSO2 PSO3K PSO2 PSO3式中：23 / 197. 6915 1018 exp1. 5128107 exp48148 exp7355. 52. 3 10 8 expT27200RTPSO3KPPSO2PO22KP 2. 2620310 5 exp76062RT35992RT420475oC475600oC11295. 3TR 1. 9873. 基本要求? 在T X 圖上，做出平衡線；至少4 條等速率線。? 以一維擬均相平推流模型為基礎，在催化劑用量最少的前提下，求總的及各段的催化劑裝量，進出口溫度、轉化率并在T X 圖上標出

3、折線。? 程序用C、 Fortran 、 BASIC語言之一編制。4. 討論? 要求的最終轉化率從98變化到99對催化劑用量的影響；? yo2 ySO2 21，SO2進口濃度在7 9之間變化，對催化劑裝量的影響。T-X 圖繪制平衡線與等反應速率線本次大作業計算程序，使用MATlab編程實現。表 1.平衡線所需數據溫度 /K987.51889.21 835.22 797.47 767.96 743.32 721.81702.40Xe0.050.10.150.20.250.30.350.4684.4167.37 60.90 634.65 618.32 601.54 583.84564.55542.

4、480.450.50.550.60.650.70.750.80.85515.1775.21 45.21 463.63 449.41 430.53 400.84373.810.90.950.950.960.970.980.990.995使用 matlab 導出的數據作平衡線圖，如圖所示。XeXe0.81 T-X 平衡線圖圖 2. 完整范圍內的T X圖圖中 Rso2 的適宜反應范圍是420600，但在更高的溫度范圍內也是會有反應的，即使反應曲線在適宜溫度范圍以外精確度低，或者反應體系發生變化，此圖權當得到更加美觀、完整的圖，以期反映出整個T-X 圖的趨勢。計算數據如下表：表 2 等反應速率線數據R

5、=-1.9489E-5R=-1.6923E-05R=-1.4462E-05R=-1.2098E-05R=-9.8175E-06T/Xso2R=-7.5949E-06R=-5.3765E-06T/Xso2R=-3.0059E-06T/Xso2T/Xso2T/Xso2T/Xso2T/ Xso2T/Xso26200.0616200.1556200.2456200.3316200.4116200.4856200.5536200.6166100.0766100.1696100.2596100.3466100.4296100.5056100.5766100.6426000.0896000.1826000.

6、2726000.3606000.4456000.5246000.5986000.6675900.1005900.1935900.2845900.3735900.4595900.5415900.6185900.6915800.1095800.2025800.2945800.3845800.4725800.5575800.6375800.7145700.1165700.2095700.3025700.3935700.4835700.5705700.6545700.7355600.1195600.2135600.3075600.4005600.4925600.5825600.6695600.7545

7、500.1185500.2145500.3095500.4035500.4985500.5915500.6825500.7725400.1125400.2105400.3075400.4045400.5015400.5975400.6925400.7875300.1005300.2005300.3005300.4005300.5005300.6005300.7005300.8005200.0795200.1835200.2875200.3915200.4955200.5995200.7055200.8115100.0485100.1575100.2665100.3755100.4845100.

8、5945100.7065100.8195000.0035000.1195000.2345000.3505000.4675000.5845000.7035000.824490-0.0604900.0644900.1884900.3134900.4394900.5664900.6954900.827480-0.147480-0.0134800.1234800.2604800.3984800.5384800.6804800.826470-0.264470-0.1184700.0314700.1844700.3394700.4964700.6564700.820460-0.981460-0.85646

9、0-0.692460-0.494460-0.264460-0.0034600.2924600.632450450-0.555450-0.1934500.337440-0.177二 .四段反應器數據的計算及優化在工業實踐中，對于任何化學反應，要保證反應在盡量高的反應速率下進行，即意味著減小反應器的體積，減少設備投資，減少催化劑的填量，即意味著減少操作費用，更快的反應速率對應著更小的反應器體積，也意味著更大的收益。對于放熱反應，反應所放熱量抑制反應向正方向進行。隨著溫度的升高，雖然正反應速率增大，但逆反應速率增大更快，整個反應便在更低的轉化率下達到平衡，放熱反應高溫對應著低轉化率。而且隨著反應溫度

10、升高，會出現反應體系高溫下出現的問題、能耗、設備要求等其他方面，一般不在過高的溫度下進行生產。從多方面考慮，但反應放熱到一定溫度，會使物料抽出與冷源換熱降低物料溫度，以達到更快的反應速率、更大的轉化率。多級反應器的級間換熱即是出于這個目的。對于四級反應器，若使反應所需催化劑用量最少。由把 Wcat分別對各段求x 和 T的微分，使其等于0；即即下一段入口溫度點的反應速率和上一段出口的反應速率相等。假設第一段入口溫度Tin(1), 進行操作線計算。1. 操作線線斜率：已知入口溫度、組成，出口組成，求出口溫度：so2 反應所放出的熱量nso2* xso2 * (- H) =M*cp* T xso2/

11、 T=M*cp/nso2* (- H)? 操作線斜率即可求得。? 由 function t1 = TXXtoT (t0,x0,x1 ) 實現? xso2/ T=1/232.30452. 求反應器出口轉化率Xout(i) 、出口溫度Tout(i) 和所用催化劑的量Wcat（i） ：已知Tout（ i-1 ） ,Xout （ i-1 ） , 求 Tin （ i ）當滿足條件（ 1）時，反應器催化劑用量在此條件下取得極值，求微分、積分可以用matlab 中的求偏導函數diff 、 積分函數int ， 亦可用其他高精度的微積分函數，但我所使用對 （1/r） 的 Ti 偏導數計算時間長，計算一萬次可花費

12、數分鐘（本人所使用的計算機） ，對偏微分的積分耗時更長，一次可達數分鐘。因此選用梯形法求微分、積分也可以用一x 的微小偏差，獲得的y 的增量，用（ （ y+ y） -y） / x即為此處導數，求積分則用梯形法，選取足夠小的步長，算出每一步長對應的面積，作為積分值。用梯形法計算對偏微分的積分，十萬次只需數秒，具有實踐性。求出口 T 見函數 function daera . 求催化劑用量見函數function dWcat .應當注意的是，Xso2滿足條件（1）時，不可大于600對應的轉化率，如果大于 600 時的轉化率，則此段出口轉化率為600對應的轉化率，出口溫度即為600。已知入口溫度，出口轉

13、化率，已知操作線斜率，即可算出出口溫度、出口時的反應速率r. 見函數 function TXXtoT。3求反應器入口溫度已知上一段出口轉化率Xout（ i-1 ） 、出口溫度Tout（i-1） 、反應速率rout（i-1） 。下一段入口轉化率Xin（ i ） =Xout（ i-1 ） , 由條件可知 r（i）=r（i-1）. 已知Xso2、 r，求 T。已知第二段入口反應速率r（Tin（2）,Xin（2）利用 function T = RXtoT2( r,x) 或者 function T = RXtoT( r,x)對于求第二段入口溫度時，如果使用割線法/牛頓法，因為在420475時，導數值很大

14、，而且有重根的情況(因為 r(Tin(2),Xin(1)=r(Tout(1)， Xin(1) ) ) ，有時候所得不是想要的解，甚至得到不到解。有時需手動改變初值。比較麻煩在已知溫度在420到Tout(i-1) 之間情況下，考慮采用變步長依次搜索法，計算 50 次，即可達到T 的 7 位小數的精確度，且屢試不爽。見函數functionRXtoT2.function dw = dWcat (xd,tin,xin)已知入口溫度，出口轉化率，已知操作線斜率，即可算出出口溫度、出口時的反應速率r. 見函數 function TXXtoT 。給定一個第一入口溫度，求得四級反應器的數據，xso2 步長 0

15、.0001 如果使用 matalb ， (從 00.98，計算 9800次)積分計算一次需要數分鐘，如果使用梯形法計算微分積分，步長選擇0.00001 (計算近10萬次) ，計算一次四級反應器數據只需不到兩秒。梯形法雖然簡單，但計算速率高了2、 3 個數量級。給定第一段進口溫度，運行一次即可在excel 中得到四段反應器數據，例子中當第一段進口溫度設為444.65時，計算四段反應器數據如下：表 3 初設第一段進口溫度設為444.65，反應器數據XinTin/XoutTout/Wcat/kg總10444.6500.66873600.0004123.7840440.7920.66873450.57

16、00.90343505.0925600.0330.90343444.8880.96075458.2049535.9140.96075422.7150.97890426.93121181.06.對第一段入口溫度的優化對第一段入口溫度Tin ( 1)的優化 , 現在溫度范圍內得到Tin ( 1)對催化劑總量和最終轉化率的影響趨勢，再在小范圍內搜索最優值。根據以上程序，使第一段入口溫度Tin(1) 從變化，得到第一段入口溫度Tin(1) 對催化劑總量和最終轉化率的影響，得到數據列表如下表 4 不同第一段進口溫度反應器數據XinTin/XoutTout/Wcat/kg10.00000420.0000.

17、72907589.3679259.0261450.0620.72907420.5510.95177472.28521700.0130.95177420.0000.97978426.50730463.7440.97978420.0000.97979420.00227.2910.00000425.0000.71928592.0937712.7456716.9020.71928425.5490.94571478.14916480.0830.94571420.0000.97939427.82432498.6440.97939420.0000.97940420.00225.4310.00000430.0

18、000.70898594.7006528.3353457.3020.70898430.5900.93891484.00412829.1230.93891420.0000.97889429.28834076.4040.97889420.0000.97890420.00223.4410.00000435.0000.69804597.1595606.1648666.8820.69804435.6260.93125489.80210178.7630.93125425.8030.97545436.07125214.9940.97545420.0000.97935420.9067666.9810.0000

19、0440.0000.68625599.4204872.5546810.1920.68625440.6800.92250495.5628202.5630.92250432.3100.97108443.59618018.0540.97108420.0000.98011422.09815717.0310.00000445.0000.66723600.0004072.9038815.8420.66723451.2630.90180505.7545455.9130.90180445.8000.95988459.2929112.3040.95988423.7130.97838428.01120174.72

20、10.00000450.0000.64570600.0003472.4420243.6720.64570458.7190.87962513.0594115.2830.87962455.5880.94828471.5385291.7940.94828441.4090.96728445.8237364.16matlab 源程序見附錄2.1根據上表數據作圖如下，圖 4最終轉化率和催化劑總量隨著Tin(1) 變化趨勢圖由上圖可以看出，在420450范圍內，存在第一段入口溫度使得最后轉化0.98，并且在440445之間。下一步從第一段入口溫度445向4400.1 。所得結果如下表 5 轉化率達標、催化劑

21、最少用量時，反應器數據總Xin Tin/ Xout Tout/ Wcat/kg10.443.8000.67239600.4254.6120.67239448.7510.90746503.3595995.5930.90746442.5550.96290455.43410707.0140.96290420.2550.98015424.26223841.2144798.42matlab 程序見附錄2.2由上表可知滿足最終轉化率0.98 時，所有催化劑量為44798.42kg. 根據上Origin 作操作折線圖如下：4 轉化率達標、催化劑最少用量時，反應器操作線三.討論：? 1.要求的最終轉化率從98

22、變化到99對催化劑用量的影響；根據平衡線圖，平衡轉化率Xe=0.99時，對應的平衡溫度Te=400.84，而平衡轉化率隨著溫度升高而減小，所以在催化劑適用范圍420600內，最終轉化率達不到0.99。 求在420600攝氏度溫度范圍內，最大轉化率時，四段反應器數據結果如下表所示：表 6 最大轉化率時，反應器數據XinTin/XoutTout/Wcat/kg總10.00000443.9800.60616600.0002301.35240754.2420.60616470.0000.87710539.7375966.96730.87710420.0000.97632445.53822586.786

23、40.97632420.0000.98368421.8949899.133計算 matlab 源程序見附錄3.1.? 2.YO2 YSO2 21，SO2進口濃度在7 9之間變化，對催化劑裝量的影響。計算 so2 進料量在保持131kmol, 求滿足轉化率0.98 條件下催化劑用量最少表 7 催化劑用量隨著進口so2 濃度變化數據Yso2, mol%0.070.080.09Tin(1)/ 450454.99448.76Wcat/kg21626.8929988.336505.3可以看出，隨著反應器進口Yso2濃度升高。所需催化劑填裝量增大。Yso2為變量，T=450，求r表7T=450，反應速率隨

24、著so2 濃度變化數據xso20.070.080.090.6-4.72E-06-4.36E-06-4.00E-060.5-4.27E-06-3.92E-06-3.57E-060.4-3.83E-06-3.49E-06-3.16E-060.3-3.39E-06-3.07E-06-2.75E-060.2-2.95E-06-2.65E-06-2.36E-060.1-2.50E-06-2.23E-06-1.96E-06Yso2為變量，T=500，求r表7T=500，反應速率隨著so2 濃度變化數據xso20.070.080.090.6-2.05E-05-1.89E-05-1.73E-050.5-1.8

25、2E-05-1.67E-05-1.52E-050.4-1.59E-05-1.45E-05-1.31E-050.3-1.37E-05-1.24E-05-1.11E-050.2-1.15E-05-0.00001-9.18E-060.1-9.29E-060.000-7.29E-06由上表可以看出，隨著 so2 入口濃度的增大，反應速率的確降低了。符合以上結論。附錄 1 T-X 圖繪制平衡線與等反應速率線matlab 代碼clc,clear,yso2=0.08;yo2=0.09;yn2=1-yso2-yo2;epsilong=-yso2*0.5;Pso20=1.01325*yso2; Xe(26)=0

26、.0; for i=1:26if i<=20Xe(i)=-0.05+i*0.05;elseif i<26Xe(i)=0.95+(i-21)*0.01;elseXe(i)=0.995;endPso2=Pso20*(1-Xe(i)/(1+epsilong*Xe(i);Po2=(yo2-yso2*0.5*Xe(i)/yso2/(1-Xe(i) *Pso2;Pso3=Pso2*Xe(i)/(1-Xe(i);Kp=Pso3/(Pso2*Po20.5);Te(i)=11295.3/log(Kp/2.26203e-5); endTe=Te-273.15;TeXe=Te;Xe;TeXe=TeXe&

27、#39;,xlswrite('Txn.xls',TeXe,'sheet2');Xso2(1)=0;for k=1:8j=2;T(1)=530+273.15;Xso2(k)=k*0.1;Rso2=TXtoRso2(T(1),Xso2(k);TXresult(1,1)=Rso2;T=1300+273.15; for i=1:90T=T-10;a=TRXtoX(T,Rso2,0.05);if a<1TXresult(j,1)=T-273.15;TXresult(j,2)=a;fprintf('rso2(%d)=%5g ,T= %0.1f ,Xso2= %

28、0.5f ,n',k,Rso2,T-273.15,a);j=j+1;elsecontinue;endendfprintf('n');if k=1xlswrite('Txn.xls',TXresult,'sheet2','C1:D90');elseif k=2xlswrite('Txn.xls',TXresult,'sheet2','E1:F90');elseif k=3xlswrite('Txn.xls',TXresult,'sheet2',&

29、#39;G1:H90');elseif k=4 xlswrite('Txn.xls',TXresult,'sheet2','I1:J90');elseif k=5 xlswrite('Txn.xls',TXresult,'sheet2','K1:L90');elseif k=6xlswrite('Txn.xls',TXresult,'sheet2','M1:N90');elseif k=7 xlswrite('Txn.xls',

30、TXresult,'sheet2','O1:P90');elseif k=8 xlswrite('Txn.xls',TXresult,'sheet2','Q1:R90');end ;clear TXresult;end附錄 2.1 計算四段反應器各段數據及優化for k=1:8Tin=415+273.15+5*k;xin=0;dx=0.00001;fprintf(2,'t0=%0.2f,x0=%0.5f',Tin-273.15,xin);xso2=0.5*dx;WCAT=0;for i=1:4sum

31、=0;sumw=0;x600=(600+273.15-Tin(i)/232.3045128+xin(i);while sum<=0da=daera(xso2,Tin(i),xin(i);dW=dWcat(xso2,Tin(i),xin(i);sum=sum+da;sumw=sumw+dW;xso2=xso2+dx;if xso2>=x600xso2=x600;break;endendXso2(i)=xso2;Wcat(i)=sumw;WCAT=WCAT+Wcat(i);Tout(i)=TXXtoT(Tin(i),xin(i), Xso2(i); if Tout(i)>873.

32、15Tout(i)=873.15; endr(i)=TXtoRso2(Tout(i), Xso2(i);fprintf(2,'xout(%d)= %0.5f ,Tout(%d)= %0.3f?Wcat(%d)=%0.5f;n',i,Xso2(i),i,Tout(i)-273.15,i,Wcat(i);Tin(i+1)=RXtoT2(r(i),Xso2(i);if Tin(i+1)<693.15Tin(i+1)=693.15;endxin(i+1)=xso2;fprintf(2,'Tin(%d)=%0.4f? ?n',i+1,Tin(i+1)-273.15)

33、;Q2result(i,:)=i,xin(i),Tin(i)-273.15,Xso2(i),Tout(i)-273.15,Wcat(i);endif k=1xlswrite('Txn.xls',Q2result,'sheet2','S1:X4');elseif k=2xlswrite('Txn.xls',Q2result,'sheet2','S6:X9');elseif k=3 xlswrite('Txn.xls',Q2result,'sheet2','S11

34、:X14');elseif k=4 xlswrite('Txn.xls',Q2result,'sheet2','S16:X19');elseif k=5 xlswrite('Txn.xls',Q2result,'sheet2','S21:X24');elseif k=6xlswrite('Txn.xls',Q2result,'sheet2','S26:X29');elseif k=7xlswrite('Txn.xls',Q2re

35、sult,'sheet2','S31:X34')endend附錄 2.2 計算四段反應器各段數據及優化for k=1:50Tin=445+273.15-0.01*k;xin=0;dx=0.00001;fprintf(2,'t0=%0.2f,x0=%0.5f',Tin-273.15,xin);xso2=0.5*dx;WCAT=0;for i=1:4sum=0;sumw=0;x600=(600+273.15-Tin(i)/ 232.3045128+xin(i);while sum<=0da=daera(xso2,Tin(i),xin(i);dW

36、=dWcat(xso2,Tin(i),xin(i);sum=sum+da;sumw=sumw+dW;xso2=xso2+dx;if xso2>=x600xso2=x600;break;endendXso2(i)=xso2;Wcat(i)=sumw;WCAT=WCAT+Wcat(i);Tout(i)=TXXtoT(Tin(i),xin(i), Xso2(i); if Tout(i)>873.15Tout(i)=873.15; endr(i)=TXtoRso2(Tout(i), Xso2(i);fprintf(2,'k=%d xout(%d)= %0.5f ,Tout(%d)=

37、 %0.3f?Wcat(%d)=%0.5f;n',k,i,Xso2(i),i,Tout(i)-273.15,i,Wcat(i);Tin(i+1)=RXtoT2(r(i),Xso2(i);if Tin(i+1)<693.15Tin(i+1)=693.15;endxin(i+1)=xso2;fprintf(2,'Tin(%d)=%0.4f? ?n',i+1,Tin(i+1)-273.15);Q2result(i,:)=i,xin(i),Tin(i)-273.15,Xso2(i),Tout(i)-273.15,Wcat(i); endif i=4&&xs

38、o2>0.98xlswrite('Txn.xls',Q2result,'sheet2','S36:X39');break;endend附錄 3 求在420600范圍內，催化劑用量最少情況下所達到的最大轉化率Maxso2=0for k=1:50Tin=444.0+273.15-0.01*k;xin=0;dx=0.00001;xso2=0.5*dx;WCAT=0;for i=1:4sum=0;sumw=0;x600=(600+273.15-Tin(i)/231.1255145+xin(i);while sum<=0da=daera(xso

39、2,Tin(i),xin(i);dW=dWcat(xso2,Tin(i),xin(i);sum=sum+da;sumw=sumw+dW;xso2=xso2+dx;if xso2>=x600xso2=x600;break;endendXso2(i)=xso2;Wcat(i)=sumw;WCAT=WCAT+Wcat(i);Tout(i)=TXXtoT(Tin(i),xin(i), Xso2(i); if Tout(i)>873.15Tout(i)=873.15; endr(i)=TXtoRso2(Tout(i), Xso2(i);Tin(i+1)=RXtoT2(r(i),Xso2(i)

40、;if Tin(i+1)<693.15Tin(i+1)=693.15;endxin(i+1)=xso2;Q2result(i,:)=i,xin(i),Tin(i)-273.15,Xso2(i),Tout(i)-273.15,Wcat(i);endif i=4Maxso2(k+1)=xso2,WCATif Maxso2(k+1)<Maxso2(k)xlswrite('Txn.xls',Q2result,'sheet2','S36:X39');break;endendend附錄 4 計算過程中所用到的函數function rso2 = T

41、XtoRso2( tx,xso2 )%UNTITLED3 Summary of this function goes here yso2=0.08;yo2=0.09;yn2=1-yo2-yso2;epsilong=-yso2*0.5;Pso20=1.01325*yso2; R=1.987;Pso2=Pso20*(1-xso2)/(1+epsilong*xso2);Po2=(yo2-yso2*0.5*xso2)/yso2/(1-xso2)*Pso2;Pso3=Pso2*xso2/(1-xso2);Kp=Pso3/(Pso2*Po20.5);Te=11295.3/log(Kp/2.26203e-5

42、);if tx-273.15<470keff=7.6915e18*exp(-76062/R/tx);elsekeff=1.5128e7*exp(-35992/R/tx);endB=48148*exp(-7355.5/tx);K=2.3e-8*exp(27200/R/tx);Kp=2.26203e-5*exp(11295.3/t x);xi=Pso3/(Kp*Pso2*sqrt(Po2);rso2=-keff*Po2*K*Pso2/Pso3*(1-xi2)/(sqrt(B+(B-1)*Pso2/Pso3)+sqrt(K*Pso2/Pso3 )2;endfunction RXT = RXtoT( r,x)t=415+273.15;T(1)=t;T(2)=t+5;f(1)=-TXtoRso2( T(1),x )+r;f(2)=-TXtoRso2( T(2),x )+r;j=2;while abs(T(j)-T(j-1)>0.000001&&j<500T(j+1)=T(j)-f(j)*(T(j)-T(j-1)/(f(j)-f(j-1);j=j+1;f(j)=-TXtoRso2(T(j),x)+r;endif j>=499fprintf('error in funct