楊向文，陳夕松，楊 衛(wèi)，蔣 宇

(1.東南大學(xué)自動(dòng)化學(xué)院，南京 210096；2.南京富島信息工程有限公司)

連續(xù)催化重整(簡(jiǎn)稱重整)是煉化企業(yè)加工石油的主流工藝之一，主要用于高辛烷值汽油、輕芳烴和氫氣的生產(chǎn)，對(duì)企業(yè)整體經(jīng)濟(jì)效益至關(guān)重要。對(duì)重整過程進(jìn)行流程模擬，可以優(yōu)化生產(chǎn)操作，發(fā)現(xiàn)過程薄弱環(huán)節(jié)，有助于提高企業(yè)效益。

流程模擬的效果取決于構(gòu)建模型的精準(zhǔn)性，對(duì)于重整過程的模擬目前多采用集總模型。構(gòu)建集總模型，一般先根據(jù)反應(yīng)物分子結(jié)構(gòu)和動(dòng)力學(xué)特性劃分集總，建立初始模型；然后對(duì)初始模型進(jìn)行參數(shù)估計(jì)，并用實(shí)驗(yàn)室試驗(yàn)數(shù)據(jù)或工業(yè)生產(chǎn)實(shí)際數(shù)據(jù)修正。雖然模型參數(shù)估計(jì)有理論方法為依據(jù)，但是理論參數(shù)估計(jì)得到的模型與工業(yè)實(shí)際生產(chǎn)過程往往存在偏差，因此需要使用實(shí)際生產(chǎn)數(shù)據(jù)對(duì)集總模型參數(shù)進(jìn)行修正。

重整反應(yīng)體系極其復(fù)雜，有近300種組分，重整過程中發(fā)生的反應(yīng)更是不計(jì)其數(shù)[1]。各反應(yīng)的速率受到反應(yīng)物濃度、反應(yīng)壓力、反應(yīng)溫度、催化劑活性等諸多因素的影響，反應(yīng)物間的相互轉(zhuǎn)化存在較強(qiáng)耦合性，不同反應(yīng)的反應(yīng)速率也存在較大差異，致使依靠傳統(tǒng)算法的梯度求解難度很大。以30個(gè)集總為例，若采用傳統(tǒng)算法進(jìn)行一次參數(shù)估計(jì)得到滿意結(jié)果，往往需要調(diào)用近萬次集總模型，參數(shù)估計(jì)效率很低。因此，優(yōu)化參數(shù)估計(jì)算法勢(shì)在必行。有學(xué)者采用變尺度法(BFGS)[2]和Marquardt法[3]等對(duì)梯度運(yùn)算進(jìn)行優(yōu)化，但隨著集總劃分不斷精細(xì)化，模型愈加復(fù)雜，梯度運(yùn)算仍然難以進(jìn)行。徐斌等[4]和王建偉等[5]將進(jìn)化算法和群智能算法等智能算法用于提高參數(shù)估計(jì)速率，智能算法雖然不依賴梯度計(jì)算、速率高，但其往往只專注于算法本身性能的提升，易出現(xiàn)局部最優(yōu)解和收斂效果不佳等問題。從生產(chǎn)工藝角度分析，集總模型存在偏差主要是因?yàn)槟Ｐ痛嬖趧傂訹6]。Dotto等[7]優(yōu)化了剛性微分方程的數(shù)值求解算法，采用Gear法對(duì)剛性集總模型進(jìn)行計(jì)算，但Gear法需要迭代多次求解，求解效率有待提升。

上述方法均未能徹底解決集總模型的剛性問題和集總模型參數(shù)估計(jì)效率較低問題，本研究通過深入分析重整生產(chǎn)工藝，采用準(zhǔn)穩(wěn)態(tài)假設(shè)簡(jiǎn)化集總模型，在遵循反應(yīng)真實(shí)變化規(guī)律的前提下去除集總模型的剛性，在保證模型計(jì)算精度的條件下有效提高集總模型參數(shù)估計(jì)效率，期望高效精準(zhǔn)地對(duì)重整過程進(jìn)行流程模擬。

1 集總模型參數(shù)估計(jì)效率的優(yōu)化

1.1 重整集總模型建立

建立重整集總模型，首先按結(jié)構(gòu)相似、動(dòng)力學(xué)特性相同的原則，對(duì)重整體系中的反應(yīng)物進(jìn)行分類，即劃分集總[8]，然后將集總作為虛擬的單一組分，建立簡(jiǎn)化的集總反應(yīng)網(wǎng)絡(luò)動(dòng)力學(xué)模型。集總劃分越細(xì)，模型越精確，計(jì)算難度越大，目標(biāo)函數(shù)越難以收斂；集總劃分過粗，則會(huì)忽視反應(yīng)物間的性能差異，導(dǎo)致模型精度不足。因此，在保證模型精度的前提下，為簡(jiǎn)便計(jì)算將重整過程中的反應(yīng)物劃分為30個(gè)集總，建立30集總反應(yīng)動(dòng)力學(xué)模型。具體的集總劃分為：按照碳數(shù)為C1～C12將烷烴劃分P1～P12共12個(gè)集總；按照碳數(shù)為C6～C12將環(huán)烷烴劃分5N6(甲基環(huán)戊烷)、6N6(環(huán)己烷)和N7～N12環(huán)烷烴共8個(gè)集總；按照碳數(shù)為C6～C12將芳烴劃分A6～A12共10個(gè)集總，其中A8包括乙苯(EB)、間二甲苯(PX)、對(duì)二甲苯(MX)、鄰二甲苯(OX)4個(gè)集總。

30集總模型的反應(yīng)網(wǎng)絡(luò)如圖1所示。由圖1可知，30集總模型的反應(yīng)網(wǎng)絡(luò)包括了7種烷烴脫氫環(huán)化生成環(huán)烷烴反應(yīng)、7種環(huán)烷烴脫氫環(huán)化生成芳烴反應(yīng)、6種芳烴氫解反應(yīng)和20種烷烴加氫裂化反應(yīng)，共計(jì)40種反應(yīng)。

圖1 重整體系30集總模型的反應(yīng)網(wǎng)絡(luò)

重整過程集總模型的建立使用實(shí)驗(yàn)室試驗(yàn)數(shù)據(jù)進(jìn)行參數(shù)估計(jì)，其與工業(yè)實(shí)際生產(chǎn)情況存在偏差。為使所建模型符合工業(yè)實(shí)際情況，對(duì)模型參數(shù)估計(jì)的修正歸結(jié)為求解40個(gè)參數(shù)修正因子(λj，j=1，2，…，40)。反應(yīng)動(dòng)力學(xué)中，反應(yīng)物的濃度變化可以用其在流向不同位置上摩爾流量變化率來描述。30個(gè)集總模型中每個(gè)集總的摩爾流量變化率常微分方程組(ODE)如式(1)所示。

(1)

式(1)中：Yi為集總i的摩爾流量，kmol/h；Y為30個(gè)集總摩爾流量集合，Y=(Y1，Y2，…，Y30)；Yi0為集總i在反應(yīng)器入口的初始摩爾流量，kmol/h；x為相對(duì)于反應(yīng)器外表面的反應(yīng)物流向坐標(biāo)，m，如圖2所示，反應(yīng)物從反應(yīng)器外表面向中心做徑向流動(dòng)；A為反應(yīng)物流向截面積，m2，對(duì)于徑向反應(yīng)器為A=2πxH，其中H為反應(yīng)器高度，m；Vc為催化劑裝填量，m3；LHSV為液時(shí)空速，h-1；r+和r-分別為集總i在所有反應(yīng)中的總生成速率和總消耗速率，kmol/h。

圖2 重整反應(yīng)器中反應(yīng)物流向示意

可逆反應(yīng)和不可逆反應(yīng)的反應(yīng)速率(r)的計(jì)算方法不同。以烷烴脫氫反應(yīng)為例，可逆反應(yīng)的反應(yīng)速率計(jì)算如式(2)所示。

r=k×(YPs-YNs/keq)

(2)

以芳烴氫解反應(yīng)為例，不可逆反應(yīng)的反應(yīng)速率計(jì)算如式(3)所示。

r=k×YAs

(3)

式(2)和式(3)中：s為碳原子數(shù)；YPs為反應(yīng)物烷烴Ps的摩爾流量，kmol/h；YNs為生成物環(huán)烷烴Ns的摩爾流量，kmol/h；YAs為反應(yīng)物芳烴As的摩爾流量，kmol/h；k為反應(yīng)速率常數(shù)，h-1，keq為可逆反應(yīng)平衡常數(shù)，k和keq的計(jì)算如式(4)～式(6)所示。

(4)

keqj=e-ΔGmol/RTj；j=1，2，…，40

(5)

(6)

1.2 準(zhǔn)穩(wěn)態(tài)假設(shè)

1.2.1集總模型中的剛性組分

進(jìn)行準(zhǔn)穩(wěn)態(tài)假設(shè)時(shí)，需要先找到模型中的剛性組分，然后再進(jìn)行穩(wěn)態(tài)假設(shè)。對(duì)于某化學(xué)反應(yīng)體系，剛性組分一般滿足式(7)描述特征，即剛性組分為體系中摩爾流量很小但流量變化率很大的組分，一般是不重要的中間產(chǎn)物。

(7)

式中，δ為較小閾值，r+-r-≈0。

1.2.2準(zhǔn)穩(wěn)態(tài)假設(shè)

在化學(xué)反應(yīng)過程中，剛性組分的摩爾流量在反應(yīng)初期變化快，反應(yīng)后期變化逐漸緩和，直至達(dá)到平衡狀態(tài)[11]；穩(wěn)態(tài)組分的摩爾流量在整個(gè)反應(yīng)過程始終處于動(dòng)態(tài)平衡狀態(tài)，不會(huì)發(fā)生較大變化。準(zhǔn)穩(wěn)態(tài)假設(shè)，即將剛性組分近似看作穩(wěn)態(tài)組分，從而簡(jiǎn)化計(jì)算、去除模型的剛性。

將剛性組分近似作為穩(wěn)態(tài)組分后，其摩爾流量基本不變，即r+-r-=0，則：

(8)

(9)

然后將式(9)代入到集總模型式(1)中，求解其濃度。集總模型從式(1)簡(jiǎn)化為式(10)，將剛性組分從ODE求解體系中剝離出來，以此降低模型的剛性和求解計(jì)算規(guī)模。

(10)

將準(zhǔn)穩(wěn)態(tài)假設(shè)簡(jiǎn)化的30集總模型與帶剛性組分的30集總模型進(jìn)行對(duì)比。首先對(duì)比使用ODE求解器獲得的兩模型中集總摩爾流量的變化曲線，分析準(zhǔn)穩(wěn)態(tài)假設(shè)模擬集總組分化學(xué)動(dòng)力學(xué)特性的準(zhǔn)確性；進(jìn)而對(duì)比兩個(gè)模型的參數(shù)估計(jì)效率，分析準(zhǔn)穩(wěn)態(tài)假設(shè)對(duì)集總模型參數(shù)估計(jì)效率的影響。

2 模型驗(yàn)證與分析

2.1 實(shí)例驗(yàn)證

利用某煉化企業(yè)重整裝置的生產(chǎn)數(shù)據(jù)對(duì)集總模型進(jìn)行實(shí)例驗(yàn)證。該重整裝置有4個(gè)反應(yīng)器(M=4)，對(duì)于每個(gè)反應(yīng)器而言，集總模型結(jié)構(gòu)相同，而模型初始值和反應(yīng)器參數(shù)不同。為避免重復(fù)，以第一反應(yīng)器為例，首先根據(jù)式(1)建立30集總模型ODE，然后使用MATLAB中ODE23S剛性求解器對(duì)集總模型ODE進(jìn)行求解，進(jìn)而用內(nèi)點(diǎn)法進(jìn)行模型參數(shù)估計(jì)，獲得第一反應(yīng)器集總模型中30個(gè)集總的摩爾流量變化情況。結(jié)果發(fā)現(xiàn)，N11和N12集總的摩爾流量變化率明顯大于其他集總，與N11和N12有關(guān)的化學(xué)反應(yīng)式如式(11)所示。

(11)

從式(11)可以看出，N11為P11轉(zhuǎn)化為A11的中間產(chǎn)物，N12為P12轉(zhuǎn)化為A12的中間產(chǎn)物。集總模型中與N11和N12有關(guān)組分的摩爾流量變化曲線如圖3所示。由圖3可見：反應(yīng)初期，相較于P11和A11，反應(yīng)物流向坐標(biāo)小于0.05 m時(shí)，N11摩爾流量急劇增大，其變化曲線的斜率極大，但是其摩爾流量比P11和A11的摩爾流量小1～2個(gè)數(shù)量級(jí)；而反應(yīng)后期(流向坐標(biāo)大于0.2 m)，N11和N12摩爾流量變化曲線的斜率逐漸減小，摩爾流量趨于穩(wěn)定。N11和N12摩爾流量的變化情況符合剛性組分的變化特征，即N11和N12為集總模型中的剛性組分。因此，將N11和N12近似看作穩(wěn)態(tài)組分，進(jìn)行準(zhǔn)穩(wěn)態(tài)假設(shè)簡(jiǎn)化。

圖3 P11，N11，A11，P12，N12，A12組分的摩爾流量變化曲線

根據(jù)式(1)和式(2)，集總組分N11和N12的凈生成速率分別由式(12)和式(13)計(jì)算。

(12)

(13)

(14)

(15)

將式(14)和式(15)代入30集總簡(jiǎn)化模型式(10)中求解。這樣就不再需要求解剛性組分N11和N12的摩爾流量，從而降低了模型的剛性，進(jìn)而采用非剛性求解器ODE23對(duì)集總模型進(jìn)行更快地求解和參數(shù)估計(jì)。準(zhǔn)穩(wěn)態(tài)假設(shè)前后P11，N11，A11，P12，N12，A12摩爾流量的變化曲線如圖4所示。

圖4 準(zhǔn)穩(wěn)態(tài)假設(shè)前后P11，N11，A11，P12，N12，A12的摩爾流量變化曲線對(duì)比●—原始曲線; ▲—準(zhǔn)穩(wěn)態(tài)假設(shè)曲線

由圖4可見：采用準(zhǔn)穩(wěn)態(tài)假設(shè)后，P11和P12摩爾流量的變化曲線與準(zhǔn)穩(wěn)態(tài)假設(shè)前一致；而N11，N12，A11，A12摩爾流量的變化趨勢(shì)不變，但與準(zhǔn)穩(wěn)態(tài)假設(shè)前存在一定偏差。這表明，采用準(zhǔn)穩(wěn)態(tài)假設(shè)較好地保持了原始集總組分的化學(xué)動(dòng)力學(xué)特性，符合反應(yīng)中組分摩爾流量的真實(shí)變化情況。準(zhǔn)穩(wěn)態(tài)假設(shè)后組分摩爾流量變化曲線中的偏差被稱為QSSA瞬時(shí)誤差，該誤差代表了計(jì)算準(zhǔn)穩(wěn)態(tài)組分摩爾流量時(shí)存在的連續(xù)擾動(dòng)[12]。反應(yīng)過程中，P11和P12先轉(zhuǎn)化為N11和N12再轉(zhuǎn)化為A11和A12。P11和P12在這6個(gè)集總中為原料集總，因此P11和P12流量計(jì)算并未受到準(zhǔn)穩(wěn)態(tài)假設(shè)后N11和N12摩爾流量計(jì)算偏差的影響；而A11和A12在這6個(gè)集總中為產(chǎn)物集總，其摩爾流量的計(jì)算直接受到N11和N12摩爾流量計(jì)算偏差的影響，導(dǎo)致其變化曲線存在一定偏差。

2.2 模型參數(shù)估計(jì)

采用準(zhǔn)穩(wěn)態(tài)假設(shè)去除模型剛性組分后，對(duì)簡(jiǎn)化模型進(jìn)行參數(shù)估計(jì)，并與含剛性組分的30集總模型的參數(shù)估計(jì)效率進(jìn)行對(duì)比。模型參數(shù)估計(jì)時(shí)，設(shè)定決策變量為40個(gè)模型參數(shù)修正因子(λi)，設(shè)定目標(biāo)函數(shù)為理論模型中的摩爾流量、溫度等參數(shù)與實(shí)際情況差值的加權(quán)總和。目標(biāo)函數(shù)的最小值為：

(16)

準(zhǔn)穩(wěn)態(tài)假設(shè)前后集總模型參數(shù)估計(jì)的目標(biāo)函數(shù)變化曲線如圖5所示，為更清晰地展示目標(biāo)函數(shù)的變化情況，對(duì)目標(biāo)函數(shù)取對(duì)數(shù)表征。準(zhǔn)穩(wěn)態(tài)假設(shè)前后集總模型參數(shù)估計(jì)的指標(biāo)如表1所示。

圖5 準(zhǔn)穩(wěn)態(tài)假設(shè)前后模型參數(shù)估計(jì)的目標(biāo)函數(shù)變化 —原始模型； —準(zhǔn)穩(wěn)態(tài)假設(shè)簡(jiǎn)化模型

表1 準(zhǔn)穩(wěn)態(tài)假設(shè)前后集總模型參數(shù)估計(jì)的指標(biāo)

從圖5和表1可見：與帶剛性的30集總模型相比，采用準(zhǔn)穩(wěn)態(tài)假設(shè)簡(jiǎn)化去除了集總模型的剛性，模型參數(shù)估計(jì)的計(jì)算速率提升了約10倍，估計(jì)效率得到顯著提高；而且，由于去除了模型的剛性組分，減小了集總組分之間摩爾流量變化率的差異，目標(biāo)函數(shù)迭代計(jì)算的最終值大幅減小，說明簡(jiǎn)化模型計(jì)算結(jié)果與實(shí)際情況的一致性比剛性模型大幅提高；同時(shí)，在基本不損失模型計(jì)算精度的情況下，簡(jiǎn)化集總模型更易擬合，目標(biāo)函數(shù)的收斂更加迅速，有效提升了模型參數(shù)估計(jì)的效率。

綜上所述，采用準(zhǔn)穩(wěn)態(tài)假設(shè)后，雖然剛性組分及部分非剛性產(chǎn)物組分的摩爾流量計(jì)算出現(xiàn)一定偏差，但各組分的摩爾流量仍保持了原來的變化趨勢(shì)，符合實(shí)際情況；而且，采用準(zhǔn)穩(wěn)態(tài)假設(shè)后模型參數(shù)估計(jì)的效率得到大幅提升，較小的計(jì)算精度損失是可接受的。

3 結(jié) 論

在連續(xù)重整過程模擬中，采用準(zhǔn)穩(wěn)態(tài)假設(shè)簡(jiǎn)化剛性集總模型可行。通過采用準(zhǔn)穩(wěn)態(tài)假設(shè)去除集總模型的剛性后，對(duì)簡(jiǎn)化集總模型進(jìn)行參數(shù)估計(jì)，雖然計(jì)算精度有所損失，但不會(huì)影響集總模型的化學(xué)動(dòng)力學(xué)特性，反應(yīng)物流量的變化曲線基本符合實(shí)際情況。

案例分析結(jié)果表明，采用準(zhǔn)穩(wěn)態(tài)假設(shè)去除了集總模型的剛性，在基本不損失模型計(jì)算精度的情況下，簡(jiǎn)化了集總模型計(jì)算，顯著提升了集總模型參數(shù)估計(jì)的效率，而且簡(jiǎn)化模型計(jì)算結(jié)果與實(shí)際情況的一致性比剛性模型大幅提高。