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

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

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

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

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

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

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

1.1 重整集總模型建立

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

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

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

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

(1)

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

圖2 重整反應器中反應物流向示意

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

r=k×(YPs-YNs/keq)

(2)

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

r=k×YAs

(3)

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

(4)

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

(5)

(6)

1.2 準穩(wěn)態(tài)假設

1.2.1集總模型中的剛性組分

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

(7)

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

1.2.2準穩(wěn)態(tài)假設

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

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

(8)

(9)

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

(10)

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

2 模型驗證與分析

2.1 實例驗證

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

(11)

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

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

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

(12)

(13)

(14)

(15)

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

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

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

2.2 模型參數(shù)估計

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

(16)

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

圖5 準穩(wěn)態(tài)假設前后模型參數(shù)估計的目標函數(shù)變化 —原始模型； —準穩(wěn)態(tài)假設簡化模型

表1 準穩(wěn)態(tài)假設前后集總模型參數(shù)估計的指標

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

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

3 結 論

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

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