1、化學(xué)反應(yīng)器理論大作業(yè)化學(xué)反應(yīng)器理論大作業(yè)二氧化硫轉(zhuǎn)換器最優(yōu)化班 級: 化研1612 學(xué) 號: 指導(dǎo)教師: 文利雄 姓 名: 閆曉宇 18 / 19二氧化硫轉(zhuǎn)換器最優(yōu)化反應(yīng)方程式：SO21/2O2=SO3 （放熱反應(yīng)）四段絕熱反應(yīng)器，級間間接換熱，常壓下反應(yīng)。1. 基礎(chǔ)數(shù)據(jù) 混合物恒壓熱容Cp 0.2549 kcal/kgK H =23135 kcal/kmol 催化劑堆密度 b554 kg/m3 進(jìn)口SO2濃度8.0 mol%，O2濃度9.0 mol%，其余為氮?dú)?處理量131 kmol SO2/hr，要求最終轉(zhuǎn)化率98 2.動力學(xué)方程式中： 3. 基本要求 在TX 圖上，做出平衡線；至少4條

2、等速率線。 以一維擬均相平推流模型為基礎(chǔ)，在催化劑用量最少的前提下，求總的及各段的催化劑裝量，進(jìn)出口溫度、轉(zhuǎn)化率并在TX 圖上標(biāo)出折線。 程序用C、Fortran、BASIC語言之一編制。4.討論 要求的最終轉(zhuǎn)化率從98變化到99對催化劑用量的影響； yo2ySO221，SO2進(jìn)口濃度在79之間變化，對催化劑裝量的影響。一T-X圖繪制平衡線與等反應(yīng)速率線 本次大作業(yè)計(jì)算程序，使用MATlab編程實(shí)現(xiàn)。 表1.平衡線所需數(shù)據(jù)溫度/K987.51 889.21 835.22 797.47 767.96 743.32 721.81 702.40 Xe0.050.10.150.20.250.30.35

3、0.4684.41 667.37 650.90 634.65 618.32 601.54 583.84 564.55 542.48 0.450.50.550.60.650.70.750.80.85515.17 475.21 475.21 463.63 449.41 430.53 400.84 373.81 0.90.950.950.960.970.980.990.995使用matlab導(dǎo)出的數(shù)據(jù)作平衡線圖，如圖所示。圖1 T-X平衡線圖圖2. 完整范圍內(nèi)的TX圖圖中Rso2的適宜反應(yīng)范圍是420600，但在更高的溫度范圍內(nèi)也是會有反應(yīng)的，即使反應(yīng)曲線在適宜溫度范圍以外精確度低，或者反應(yīng)體系發(fā)生

4、變化，此圖權(quán)當(dāng)?shù)玫礁用烙^、完整的圖，以期反映出整個(gè)T-X圖的趨勢。計(jì)算數(shù)據(jù)如下表：表2 等反應(yīng)速率線數(shù)據(jù)R=-1.9489E-5R=-1.6923E-05R=-1.4462E-05R=-1.2098E-05R=-9.8175E-06R=-7.5949E-06R=-5.3765E-06R=-3.0059E-06T/ Xso2T/ Xso2T/ Xso2T/ Xso2T/ Xso2T/ Xso2T/ Xso2T/ Xso26200.0616200.1556200.2456200.3316200.411 6200.485 6200.553 6200.6166100.0766100.1696100.

5、2596100.3466100.429 6100.505 6100.576 6100.6426000.0896000.1826000.2726000.3606000.445 6000.524 6000.598 6000.6675900.1005900.1935900.2845900.3735900.459 5900.541 5900.618 5900.6915800.1095800.2025800.2945800.3845800.472 5800.557 5800.637 5800.7145700.1165700.2095700.3025700.3935700.483 5700.570 570

6、0.654 5700.7355600.1195600.2135600.3075600.4005600.492 5600.582 5600.669 5600.7545500.1185500.2145500.3095500.4035500.498 5500.591 5500.682 5500.7725400.1125400.2105400.3075400.4045400.501 5400.597 5400.692 5400.7875300.1005300.2005300.3005300.4005300.500 5300.600 5300.700 5300.8005200.0795200.18352

7、00.2875200.3915200.495 5200.599 5200.705 5200.8115100.0485100.1575100.2665100.3755100.484 5100.594 5100.706 5100.8195000.0035000.1195000.2345000.3505000.467 5000.584 5000.703 5000.824490-0.0604900.0644900.1884900.3134900.439 4900.566 4900.695 4900.827480-0.147480-0.0134800.1234800.2604800.398 4800.5

8、38 4800.680 4800.826470-0.264470-0.1184700.0314700.1844700.339 4700.496 4700.656 4700.820460-0.981460-0.856460-0.692460-0.494460-0.264460-0.0034600.292 4600.632450450-0.555450-0.1934500.337440-0.177二.四段反應(yīng)器數(shù)據(jù)的計(jì)算及優(yōu)化在工業(yè)實(shí)踐中，對于任何化學(xué)反應(yīng)，要保證反應(yīng)在盡量高的反應(yīng)速率下進(jìn)行，即意味著減小反應(yīng)器的體積，減少設(shè)備投資，減少催化劑的填量，即意味著減少操作費(fèi)用，更快的反應(yīng)速率對應(yīng)著更小的

9、反應(yīng)器體積，也意味著更大的收益。對于放熱反應(yīng)，反應(yīng)所放熱量抑制反應(yīng)向正方向進(jìn)行。隨著溫度的升高，雖然正反應(yīng)速率增大，但逆反應(yīng)速率增大更快，整個(gè)反應(yīng)便在更低的轉(zhuǎn)化率下達(dá)到平衡，放熱反應(yīng)高溫對應(yīng)著低轉(zhuǎn)化率。而且隨著反應(yīng)溫度升高，會出現(xiàn)反應(yīng)體系高溫下出現(xiàn)的問題、能耗、設(shè)備要求等其他方面，一般不在過高的溫度下進(jìn)行生產(chǎn)。從多方面考慮，但反應(yīng)放熱到一定溫度，會使物料抽出與冷源換熱降低物料溫度，以達(dá)到更快的反應(yīng)速率、更大的轉(zhuǎn)化率。多級反應(yīng)器的級間換熱即是出于這個(gè)目的。對于四級反應(yīng)器，若使反應(yīng)所需催化劑用量最少。由把Wcat分別對各段求x和T的微分，使其等于0；即即下一段入口溫度點(diǎn)的反應(yīng)速率和上一段出口的反應(yīng)

10、速率相等。假設(shè)第一段入口溫度Tin(1),進(jìn)行操作線計(jì)算。1. 操作線線斜率：已知入口溫度、組成，出口組成，求出口溫度：so2反應(yīng)所放出的熱量nso2*xso2 * (-H) =M*cp*Txso2/ T= M*cp/ nso2* (-H) 操作線斜率即可求得。 由function t1 = TXXtoT (t0,x0,x1 )實(shí)現(xiàn) xso2/ T=1/232.30452.求反應(yīng)器出口轉(zhuǎn)化率Xout(i)、出口溫度Tout(i)和所用催化劑的量Wcat(i)：已知Tout（i-1）,Xout（i-1）,求Tin（i）當(dāng)滿足條件 （1）時(shí)，反應(yīng)器催化劑用量在此條件下取得極值，求微分、積分可以用m

11、atlab中的求偏導(dǎo)函數(shù)diff、積分函數(shù)int，亦可用其他高精度的微積分函數(shù)，但我所使用對(1/r)的Ti偏導(dǎo)數(shù)計(jì)算時(shí)間長，計(jì)算一萬次可花費(fèi)數(shù)分鐘（本人所使用的計(jì)算機(jī)），對偏微分的積分耗時(shí)更長，一次可達(dá)數(shù)分鐘。因此選用梯形法求微分、積分也可以用一x的微小偏差，獲得的y的增量，用(（y+y）-y）/x即為此處導(dǎo)數(shù)，求積分則用梯形法，選取足夠小的步長，算出每一步長對應(yīng)的面積，作為積分值。用梯形法計(jì)算對偏微分的積分，十萬次只需數(shù)秒，具有實(shí)踐性。求出口T見函數(shù)function daera.求催化劑用量見函數(shù) function dWcat.應(yīng)當(dāng)注意的是，Xso2滿足條件（1）時(shí)，不可大于600對應(yīng)的轉(zhuǎn)

12、化率，如果大于600時(shí)的轉(zhuǎn)化率，則此段出口轉(zhuǎn)化率為600對應(yīng)的轉(zhuǎn)化率，出口溫度即為600。已知入口溫度，出口轉(zhuǎn)化率，已知操作線斜率，即可算出出口溫度、出口時(shí)的反應(yīng)速率r.見函數(shù) function TXXtoT。 3求反應(yīng)器入口溫度已知上一段出口轉(zhuǎn)化率Xout（i-1）、出口溫度Tout(i-1)、反應(yīng)速率rout(i-1)。下一段入口轉(zhuǎn)化率Xin（i）=Xout（i-1）,由條件可知r(i)=r(i-1). 已知Xso2、r，求T。已知第二段入口反應(yīng)速率r(Tin(2),Xin(2) 利用function T = RXtoT2( r,x)（變步長搜索） 或者 function T = RXto

13、T( r,x)（割線法）對于求第二段入口溫度時(shí)，如果使用割線法/牛頓法，因?yàn)樵?20475時(shí)，導(dǎo)數(shù)值很大，而且有重根的情況（因?yàn)?r(Tin(2),Xin(1)=r(Tout(1)，Xin(1)），有時(shí)候所得不是想要的解，甚至得到不到解。有時(shí)需手動改變初值。比較麻煩在已知溫度在420到Tout(i-1)之間情況下，考慮采用變步長依次搜索法，計(jì)算50次，即可達(dá)到T的7位小數(shù)的精確度，且屢試不爽。見函數(shù) function RXtoT2.function dw = dWcat (xd,tin,xin) 已知入口溫度，出口轉(zhuǎn)化率，已知操作線斜率，即可算出出口溫度、出口時(shí)的反應(yīng)速率r.見函數(shù) funct

14、ion TXXtoT。給定一個(gè)第一入口溫度，求得四級反應(yīng)器的數(shù)據(jù)，xso2步長0.0001如果使用matalb， （從00.98，計(jì)算9800次）積分計(jì)算一次需要數(shù)分鐘，如果使用梯形法計(jì)算微分積分，步長選擇0.00001（計(jì)算近10萬次），計(jì)算一次四級反應(yīng)器數(shù)據(jù)只需不到兩秒。梯形法雖然簡單，但計(jì)算速率高了2、3個(gè)數(shù)量級。給定第一段進(jìn)口溫度，運(yùn)行一次即可在excel中得到四段反應(yīng)器數(shù)據(jù)，例子中當(dāng)?shù)谝欢芜M(jìn)口溫度設(shè)為444.65時(shí)，計(jì)算四段反應(yīng)器數(shù)據(jù)如下：表3初設(shè)第一段進(jìn)口溫度設(shè)為444.65，反應(yīng)器數(shù)據(jù)XinTin/XoutTout/Wcat/kg總10444.6500.66873600.0004

15、123.7840440.7920.66873450.5700.90343505.0925600.0330.90343444.8880.96075458.2049535.9140.96075422.7150.97890426.93121181.06.對第一段入口溫度的優(yōu)化對第一段入口溫度Tin（1）的優(yōu)化,現(xiàn)在溫度范圍內(nèi)得到Tin（1）對催化劑總量和最終轉(zhuǎn)化率的影響趨勢，再在小范圍內(nèi)搜索最優(yōu)值。根據(jù)以上程序，使第一段入口溫度Tin(1)從變化，得到第一段入口溫度Tin(1)對催化劑總量和最終轉(zhuǎn)化率的影響，得到數(shù)據(jù)列表如下表4不同第一段進(jìn)口溫度反應(yīng)器數(shù)據(jù)XinTin/XoutTout/Wcat/k

16、g10.00000420.0000.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.002

17、25.4310.00000430.0000.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

18、.9067666.9810.00000440.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.97

19、838428.01120174.7210.00000450.0000.64570600.0003472.4420243.6720.64570458.7190.87962513.0594115.2830.87962455.5880.94828471.5385291.7940.94828441.4090.96728445.8237364.16matlab源程序見附錄2.1。根據(jù)上表數(shù)據(jù)作圖如下，最優(yōu)值區(qū)間圖4最終轉(zhuǎn)化率和催化劑總量隨著Tin(1)變化趨勢圖由上圖可以看出，在420450范圍內(nèi)，存在第一段入口溫度使得最后轉(zhuǎn)化率大于0.98，并且在440445之間。下一步從第一段入口溫度445向440

20、依次搜索，步長為0.1。所得結(jié)果如下表5轉(zhuǎn)化率達(dá)標(biāo)、催化劑最少用量時(shí)，反應(yīng)器數(shù)據(jù)XinTin/XoutTout/Wcat/kg總10. 443.800 0.67239 600. 4254.61 44798.4220.67239 448.751 0.90746 503.359 5995.59 30.90746 442.555 0.96290 455.434 10707.01 40.96290 420.255 0.98015 424.262 23841.21 matlab程序見附錄2.2由上表可知滿足最終轉(zhuǎn)化率0.98時(shí)，所有催化劑量為44798.42kg.根據(jù)上表數(shù)據(jù)用Origin作操作折線圖如

21、下：圖4轉(zhuǎn)化率達(dá)標(biāo)、催化劑最少用量時(shí)，反應(yīng)器操作線三.討論： 1.要求的最終轉(zhuǎn)化率從98變化到99對催化劑用量的影響；根據(jù)平衡線圖，平衡轉(zhuǎn)化率Xe=0.99時(shí)，對應(yīng)的平衡溫度Te=400.84，而平衡轉(zhuǎn)化率隨著溫度升高而減小，所以在催化劑適用范圍420600內(nèi)，最終轉(zhuǎn)化率達(dá)不到0.99。求在420600攝氏度溫度范圍內(nèi)，最大轉(zhuǎn)化率時(shí)，四段反應(yīng)器數(shù)據(jù)結(jié)果如下表所示：表6最大轉(zhuǎn)化率時(shí)，反應(yīng)器數(shù)據(jù)XinTin/XoutTout/Wcat/kg總10.00000 443.980 0.60616 600.000 2301.352 40754.2420.60616 470.000 0.87710 539.

22、737 5966.967 30.87710 420.000 0.97632 445.538 22586.786 40.97632 420.000 0.98368 421.894 9899.133 計(jì)算matlab源程序見附錄3.1. 2.YO2YSO221，SO2進(jìn)口濃度在79之間變化，對催化劑裝量的影響。計(jì)算so2進(jìn)料量在保持131kmol,求滿足轉(zhuǎn)化率0.98條件下催化劑用量最少表7催化劑用量隨著進(jìn)口so2濃度變化數(shù)據(jù)Yso2, mol%0.070.080.09Tin(1)/450454.99448.76Wcat/kg21626.8929988.336505.3可以看出，隨著反應(yīng)器進(jìn)口Ys

23、o2濃度升高。所需催化劑填裝量增大。以Yso2為變量，T=450，求r表7 T=450，反應(yīng)速率隨著so2濃度變化數(shù)據(jù)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-06以Yso2為變量，T=500，求r表7 T=500，反應(yīng)速率隨著so2濃度變化數(shù)據(jù)xso2

24、0.070.080.090.6-2.05E-05-1.89E-05-1.73E-050.5-1.82E-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入口濃度的增大，反應(yīng)速率的確降低了。符合以上結(jié)論。附錄1 T-X圖繪制平衡線與等反應(yīng)速率線matlab代碼clc,clear,yso2=0.08;yo2=0.09;yn2=1-yso2-yo2;eps

25、ilong=-yso2*0.5;Pso20=1.01325*yso2; Xe(26)=0.0;for i=1:26 if i=20 Xe(i)=-0.05+i*0.05;elseif i26 Xe(i)=0.95+(i-21)*0.01;else Xe(i)=0.995; end Pso2=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-

26、5);endTe=Te-273.15;TeXe=Te;Xe;TeXe=TeXe,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:90 T=T-10; a=TRXtoX(T,Rso2,0.05); if a1 TXresult(j,1)=T-273.15;TXresult(j,2)=a; fprintf(rso2(%d)=%5g ,T=

27、%0.1f ,Xso2= %0.5f ,n,k,Rso2,T-273.15,a); j=j+1; else continue; endendfprintf(n);if k=1 xlswrite(Txn.xls,TXresult,sheet2,C1:D90);elseif k=2 xlswrite(Txn.xls,TXresult,sheet2,E1:F90);elseif k=3 xlswrite(Txn.xls,TXresult,sheet2,G1:H90);elseif k=4 xlswrite(Txn.xls,TXresult,sheet2,I1:J90);elseif k=5 xlsw

28、rite(Txn.xls,TXresult,sheet2,K1:L90);elseif k=6 xlswrite(Txn.xls,TXresult,sheet2,M1:N90);elseif k=7 xlswrite(Txn.xls,TXresult,sheet2,O1:P90);elseif k=8 xlswrite(Txn.xls,TXresult,sheet2,Q1:R90);end ;clear TXresult;end附錄2.1 計(jì)算四段反應(yīng)器各段數(shù)據(jù)及優(yōu)化for k=1:8Tin=415+273.15+5*k;xin=0;dx=0.00001;fprintf(2,t0=%0.2f,

29、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.+xin(i);while sum=x600 xso2=x600; break; endend Xso2(i)=xso2;Wcat(i)=sumw;WCAT=WCAT+Wcat(i); Tout(i)=TXXtoT(Tin(i),xin(i), Xso2(i); if Tout(i)873.15 Tout(i)=873.15; end r(i)=TXtoRso2(Tout(i), Xso2(i); fpri

30、ntf(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.15 Tin(i+1)=693.15; end xin(i+1)=xso2; fprintf(2,Tin(%d)=%0.4fn,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 k=

31、1 xlswrite(Txn.xls,Q2result,sheet2,S1:X4);elseif k=2 xlswrite(Txn.xls,Q2result,sheet2,S6:X9);elseif k=3 xlswrite(Txn.xls,Q2result,sheet2,S11:X14);elseif k=4 xlswrite(Txn.xls,Q2result,sheet2,S16:X19);elseif k=5 xlswrite(Txn.xls,Q2result,sheet2,S21:X24);elseif k=6 xlswrite(Txn.xls,Q2result,sheet2,S26:

32、X29);elseif k=7 xlswrite(Txn.xls,Q2result,sheet2,S31:X34)endend附錄2.2 計(jì)算四段反應(yīng)器各段數(shù)據(jù)及優(yōu)化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.+xin(i);while sum=x600 xso2=x600; break; endend Xso2(i

33、)=xso2;Wcat(i)=sumw;WCAT=WCAT+Wcat(i); Tout(i)=TXXtoT(Tin(i),xin(i), Xso2(i); if Tout(i)873.15 Tout(i)=873.15; end r(i)=TXtoRso2(Tout(i), Xso2(i); fprintf(2,k=%d xout(%d)= %0.5f ,Tout(%d)= %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)0.98 xl

34、swrite(Txn.xls,Q2result,sheet2,S36:X39); break;endend附錄3 求在420600范圍內(nèi)，催化劑用量最少情況下所達(dá)到的最大轉(zhuǎn)化率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.+xin(i);while sum=x600 xso2=x600; break; endend Xso2(i)=xso2;Wcat(i)=sumw;WCAT=WCAT+

35、Wcat(i); Tout(i)=TXXtoT(Tin(i),xin(i), Xso2(i); if Tout(i)873.15 Tout(i)=873.15; end r(i)=TXtoRso2(Tout(i), Xso2(i); Tin(i+1)=RXtoT2(r(i),Xso2(i); if Tin(i+1)693.15 Tin(i+1)=693.15; end xin(i+1)=xso2; Q2result(i,:)=i,xin(i),Tin(i)-273.15,Xso2(i),Tout(i)-273.15,Wcat(i);end if i=4 Maxso2(k+1)=xso2,WCA

36、T if Maxso2(k+1)Maxso2(k) xlswrite(Txn.xls,Q2result,sheet2,S36:X39); break; endendend附錄4 計(jì)算過程中所用到的函數(shù)function rso2 = TXtoRso2( tx,xso2 )%UNTITLED3 Summary of this function goes hereyso2=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);P

37、o2=(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);if tx-273.150.&j=499 fprintf(error in function RXtoT,n) RXT=nan;else RXT=T(j);endend function Xr = TRXtoX(t3, R3,x3 )%TRTOX Summary of this function goes here x3so3(1)=x3; x3so3(2)=x3+0.001; if isnan