趙 波， 王 闊, 隋成玉

(1.遼寧石油化工大學 理學院，遼寧 撫順 113001； 2.撫順石油化工研究院，遼寧 撫順 113001)

三葉草型催化劑體相溫度分布數(shù)值模擬

趙 波1， 王 闊2, 隋成玉1

(1.遼寧石油化工大學 理學院，遼寧 撫順 113001； 2.撫順石油化工研究院，遼寧 撫順 113001)

以真實三葉草型加氫裂化催化劑三維體相環(huán)境為計算實體，以模擬工業(yè)運行溫度的函數(shù)作為邊界條件，采用無網(wǎng)格數(shù)值方法求解傅立葉傳熱方程，并通過計算結果分析加氫裂化催化反應過程中外界溫度波動對催化劑內部溫度分布的影響。結果表明，實際反應過程中催化劑內部并非等溫反應，在加氫裂化反應過程中外界的反應溫度波動以及催化劑內部的熱點分布對催化劑團簇內部的溫度場分布有一定的影響。催化劑體相內部的平均溫度也隨著反應體系放熱情況、催化劑粒徑、原料油密度、反應空速以及催化劑內部熱點分布情況的不同而有所變化。

加氫裂化； 傳熱； 無網(wǎng)格； 數(shù)值模擬； 分布計算

加氫裂化反應是一種強放熱反應，溫度和熱量又是該類反應最重要的影響因素。關于加氫裂化反應體系的能量恒算和節(jié)能研究的相關文獻眾多[1-2]，但研究一般僅限于較大空間尺度或單純反應熱層面[1-3]。在科學實驗和工業(yè)生產(chǎn)實踐中，對于小型實驗裝置而言，整個反應系統(tǒng)可以被視為等溫反應系統(tǒng)；而對工業(yè)裝置而言反應系統(tǒng)可能被視為一個比較嚴格的絕熱系統(tǒng)。在反應過程中，反應器內部的溫度分布必然會影響到催化劑床層不同位置的反應速率，進而影響到相關位置的熱量釋放。對于反應系統(tǒng)尤其是與工業(yè)裝置緊密相關的絕熱系統(tǒng)而言，反應產(chǎn)生熱量的傳遞必然引起體系溫度分布的重新調整。這一溫度的正反饋過程在加氫裂化工業(yè)生產(chǎn)中十分重要。

審視整個加氫裂化溫度正反饋過程不難發(fā)現(xiàn)，在絕熱反應體系中整個過程的熱量正是來自于催化劑顆粒自身。在整個反應過程中，每個催化劑顆粒都相當于反應體系的一個“微熱源”。這些微熱源的溫度分布情況直接或間接地影響反應體系的質量傳遞、熱量傳遞、反應進程以及產(chǎn)品分布。催化劑體相內部的溫度分布對于催化劑的壽命也有著極大的影響。因此對于催化劑內部體相溫度分布的描述意義十分重大。

雖然對于催化劑內部體相溫度分布的研究十分重要，但這一過程十分困難。傳統(tǒng)的煉化工業(yè)一般采用熱電偶進行溫度測定，但測定結果只能得到反應器內部的宏觀點溫度近似值，無法對反應過程的實時微觀熱效應進行分析。雖然近年來“非接觸式”的紅外測溫技術得到廣泛發(fā)展[4]，但由于催化劑內部微小而復雜的形態(tài)結構，目前沒有任何溫度測量設備能夠在如此微觀尺度進行溫度測量；同時，由于催化劑外觀幾何形態(tài)以及反應溫度條件的復雜性導致了常規(guī)導熱模型體系的經(jīng)驗公式難于適用。因此，基于傅立葉傳熱定律的數(shù)值模擬方法成為解決這一類問題的有效途徑。

1 理論及計算部分

1.1 傅立葉傳熱方程簡介

對于一般的傳熱過程，當存在穩(wěn)定熱源時，其溫度對于空間的分布形態(tài)可由一個泊松型的穩(wěn)態(tài)傅立葉傳熱方程確立。對于一個三維體系，其一般的數(shù)學形式如式(1)所示。

(1)

式中，F(xiàn)(x,y,z)表示體系內溫度對于空間的函數(shù)，當空間點處于體相內部v時，整個體系由一個泊松型的偏微分方程描述，其中qv(x,y,z)是體系的熱源函數(shù)，該函數(shù)表達當體系存在熱源或熱匯時，該熱源或熱匯對于空間的分布形式；λ表示計算體系的導熱系數(shù)，由該體系組成材料的物化性質決定。而當溫度點處于體相表面σ時，其溫度可以表達為第一型迪里克萊邊界條件描述，即將其表面溫度分布表達為函數(shù)G(x，y，z)。對于極少具有比較簡單的熱源或熱匯函數(shù)形式以及比較規(guī)則的幾何外觀形態(tài)和相對簡單的邊界條件的導熱體系，函數(shù)F(x,y,z)存在解析形式，同時拉普拉斯型方程的初值問題的解是不穩(wěn)定的[5]。對于絕大多數(shù)較為復雜的實際工程傳熱體系，相關模型方程體系的適定性問題討論極為復雜，因此，式(1)只能由數(shù)值方法進行求解。

1.2 加氫反應過程熱量衡算及相關方程參數(shù)的確立

在一般的加氫裂化反應過程中，催化劑體系內部的溫度分布可由三維的傅立葉傳熱方程描述，但實際工業(yè)反應條件下熱源函數(shù)以及催化劑表面溫度分布十分復雜。因此，本文對三維的傅立葉傳熱方程采用無網(wǎng)格數(shù)值方法進行求解。

一般而言，整個加氫裂化過程的精確反應熱計算比較困難，根據(jù)已有經(jīng)驗及相關文獻[6-8]，加氫裂化反應放熱范圍為200.0～500.0 kJ/kg，而原料油密度為0.85～0.95 kg/L，反應過程的體積空速為0.9～1.5 h-1。在此假定反應放熱為300.0 kJ/kg，油品密度為0.85 kg/L，則加工1 L油品放熱量為255.0 kJ/L，假定相關反應空速為0.9 h-1，則1 L催化劑1 h放熱為229.5 kJ/(L·h)，即其單位體積熱功率為2.295×10-4W/mm3。而根據(jù)相關文獻[8]，一般催化劑的導熱系數(shù)為0.27～0.32 W/(m·K)，即導熱系數(shù)為(2.7～3.2)×10-4W/(mm·K)，相關熱功率與導熱系數(shù)的比值為0.850 K/mm2。

1.3 無網(wǎng)格數(shù)值方法簡介及計算過程介紹

較為復雜的工程偏微分方程目前主要采用數(shù)值方法進行求解，常規(guī)的數(shù)值方法包括有限差分方法(FDM)、有限元方法(FEM)。

有限差分方法應用最早，但其對于復雜邊界條件的適應性很差。因此，在現(xiàn)代工程計算應用中受到了一定的限制。有限元方法自20世紀50年代開始興起，堪稱數(shù)值方法領域的一個里程碑。該方法的主導思想是將一個復雜形態(tài)的連續(xù)體劃分成有限個單元，各個單元通過一定的規(guī)則連接起來。正是基于這一指導思想，使得有限元方法在處理復雜幾何形體的各種線性和非線性問題中表現(xiàn)出很強的靈活性和魯棒性。因此其在工程領域得到廣泛應用，同時出現(xiàn)了大量的商業(yè)軟件包。

然而有限元方法也有其固有的缺點和局限性，例如有限元計算需要在計算前對計算體系進行相關的網(wǎng)格劃分，而形成相關計算網(wǎng)格的計算成本比較高。特別是對于復雜的三維體系而言，生成復雜的高質量三維網(wǎng)格相當困難[9-11]。產(chǎn)生這類問題的根源在于在有限元計算過程中使用了生成網(wǎng)格單元的相關連接信息。因此，形成一種不使用網(wǎng)格單元與網(wǎng)格劃分的方法十分重要。因此，近年來無網(wǎng)格類計算方法應運而生。

本文嘗試使用無網(wǎng)格計算方法對三維傅立葉傳熱方程進行數(shù)值求解，進而實現(xiàn)對具有不同幾何形態(tài)的圓柱型及三葉草型催化劑體相溫度空間分布形式的描述和分析。

整個計算的大略求解過程為：

首先選擇具有局域性質的局域函數(shù)，該局域函數(shù)用多項式函數(shù)為基底函數(shù)，相關基底函數(shù)可以表示為三維二次基底函數(shù)，如式(2)所示：

PT={1，x，xy，xz，y，yz，z，x2，y2，z2}

(2)

其次選擇三次樣條權重函數(shù)Wx(R)，如式(3)所示：

(3)

其中，

(4)

式中，xi為計算結點，x為采用點，Rω，x為定義的權函數(shù)支持域尺寸。

按此規(guī)則定義函數(shù)Wy(R)以及Wz(R)，最終定義權函數(shù)為：

(5)

依據(jù)已定義的基底函數(shù)以及相關權函數(shù)，分別形成兩類矩陣A及B：

(6)

(7)

定義后續(xù)無網(wǎng)格計算使用的形函數(shù)，如式(8)所示：

(8)

然后將生成的形函數(shù)分別對x、y以及z進行偏導求取，并將原函數(shù)及求導結果依次代入式(1)中的控制方程及邊界條件方程，形成所謂的“力矩陣”和“載荷矩陣”。

最后用線性最小二乘法求解力矩陣對于載荷矩陣的回歸系數(shù)，即得到計算結點函數(shù)的近似數(shù)值。再通過計算結點函數(shù)的數(shù)值重構整個采樣結點函數(shù)數(shù)值，可近似獲得原偏微分方程組體系的近似數(shù)值解。在整個計算過程中，計算節(jié)點與采樣節(jié)點的數(shù)目可以相同也可以不同，即力矩陣的行數(shù)與列數(shù)可以相同，也可以不同。這種設定可以比較靈活地適應具體的工程問題。

整個計算過程流程如圖1所示。

1.4 計算參數(shù)及邊界條件的設定

三葉草型催化劑能夠得以應用主要由于該幾何外形的催化劑相對于傳統(tǒng)的圓柱型催化劑具有較大的宏觀外表面積，以及在反應過程中可能擁有較小的系統(tǒng)壓降，但是由于三葉草型催化劑的外觀較傳統(tǒng)的圓柱型催化劑復雜，因此其在相關計算過程中的節(jié)點劃分更為復雜。三葉草型催化劑的截面外觀如圖2所示。

圖1 計算過程流程

圖2 三葉草型催化劑的截面外觀

從圖2可以看出，三葉草型催化劑包含兩個半徑，較小的半徑為三個相切的小圓柱體的半徑r，另一個半徑則是與三個小圓柱體相切的外接大圓的半徑R。相對于傳統(tǒng)的半徑R相同的圓柱型催化劑，三葉草型催化劑具有更小的體積和更大的外表面積。本文同時對半徑為2 mm、高度為10 mm的圓柱型催化劑和相同大圓半徑的三葉草型催化劑進行內部溫度場計算，并分析。

由于圓柱型催化劑幾何實體形狀的對稱性比較好，在計算過程中，可直接將直角坐標系映射為柱坐標系。分別在柱體的半徑方向取6個徑向采樣點，角向方向取21個采樣點，同時在軸向方向取41個采樣點，整個圓柱體共計生成采樣點4 346個，其中表面邊界采樣點421個。在整個計算過程中定義的權函數(shù)支持域尺寸為2，計算點及采樣點分布如圖3所示。

由于三葉草型催化劑的截面外觀相對復雜，同時對稱性遠比圓柱型催化劑差，其在赤道面角向方向具有更為豐富的信息。為保持采樣的近似一致性，對于相同高度及大圓半徑的三葉草型催化劑，采樣包含11個徑向采樣、42個赤道角向采樣以及41個軸向采樣。計算過程中包含9 840個采樣點，其中7 800個內部采樣點，2 040個邊界點以及4 920個計算點，其空間分布如圖4所示。在整個計算過程中定義的權函數(shù)支持域尺寸為2。

(a) 軸向

(b) 徑向圖3 圓柱型徑向及軸向采樣點及計算點分布

(a) 軸向

(b) 徑向圖4 三葉草型徑向及軸向采樣點及計算點分布

整個計算過程中模擬了在外界溫度為370 ℃的等溫環(huán)境下，熱點均勻分布的不同形狀催化劑的溫度分布情況。

1.5 模型所涉及的計算體系

整個計算過程使用的是分布式的高密度計算系統(tǒng)。其整套設備包含常規(guī)處理器計算機5臺，共計160核心。整個計算集群包含顯示體系和計算體系。顯示體系通過KVM切換器由遠程主機控制顯示。計算體系則由24口交換機完成計算數(shù)據(jù)的傳輸和交互。計算體系的每個計算節(jié)點可以單獨并行使用，也可以整體分布使用。該系統(tǒng)十分適合高密度計算占優(yōu)的系統(tǒng)體系。整個計算過程采用并行計算方式進行，形函數(shù)及其偏導數(shù)計算以及力矩陣及載荷矩陣的組裝過程采用多機并行計算[11-12]，共有體系的160核心并行參與。

2 結果與討論

2.1 等溫反應條件下熱點均勻圓柱型催化劑內部溫度場分析

在不同反應工況及催化劑尺度條件下，計算熱點分布均勻的圓柱型催化劑在外表面溫度為370 ℃下進行反應時的催化劑內部溫度分布情況，計算所獲取催化劑體系內部的最高溫度、平均溫度以及溫度分布的標準差結果見表1。

表1 熱點分布均勻的圓柱型催化劑表面溫度370 ℃及不同工況下的內部溫度

由表1可知，相關裂化反應的反應熱、原料油密度、反應空速以及催化劑半徑及長度均直接或間接影響催化劑內部的溫度分布。從單因素角度分析可以看出，催化劑顆粒內部最高溫度及平均溫度均隨反應熱、原料油密度、反應空速以及催化劑半徑、催化劑的軸向長度的增加而增大。就其影響效果而言，反應熱及催化劑顆粒半徑對于催化劑最高溫度及平均溫度的影響要顯著地大于催化劑長度、反應空速及原料油密度的影響。

表1所對應的相關不同工藝條件其催化劑內部的溫度分布比較相似，因此將7號計算對應的三維溫度場分布做二維軸向及徑向切面圖，結果如圖5所示。

(a) 整體

(b) 軸向

(c) 徑向圖5 序號7條件圓柱型催化劑內部溫度分布切片圖

對于熱點分布均勻的等溫反應條件，雖然不同計算工況對應的最高溫度及平均溫度數(shù)值有所不同，但其溫度分布形式基本類似?？梢哉J為其溫度分布的等勢面在徑向方向是嚴格的圓型分布，而在軸向剖面上比較近似橢圓分布，在催化劑的外表面溫度為370 ℃的催化劑內部溫度逐漸增加，其中徑向方向溫度梯度大于軸向方向。圓柱型催化劑的中心溫度最高，可以認為雖然在理論上催化劑內部并非等溫環(huán)境，但其內部的反應環(huán)境與等溫環(huán)境比較類似。但當催化劑半徑及長度較大時，其催化劑圓柱體中心最高溫度與比表面溫度相差較大，在該催化劑內部的反應相對于等溫反應偏離較大。因此，僅就熱效應而言，制備圓柱型催化劑時其半徑及長度一般不宜過大。

2.2 等溫反應條件下熱點均勻三葉草型催化劑內部溫度場分析

在不同反應工況及催化劑尺度條件下，計算熱點分布均勻的三葉草型催化劑在外表面溫度370 ℃下進行反應時的催化劑內部溫度分布情況，計算所獲取催化劑體系內部的最高溫度、平均溫度以及溫度分布的標準差結果見表2。

由表2可知，與圓柱型催化劑相同，相關裂化反應的反應熱、原料油密度、反應空速以及催化劑半徑均直接或間接地影響催化劑內部的溫度分布。從單因素角度分析可以看出，催化劑顆粒內部最高溫度及平均溫度均隨反應熱、原料油密度、反應空速以及催化劑半徑、催化劑的軸向長度的增加而增大。就其影響效果而言，反應熱及催化劑顆粒半徑對于催化劑最高溫度及平均溫度的影響顯著大于催化劑長度、反應空速及原料油密度的影響。通過數(shù)據(jù)對比可以發(fā)現(xiàn)，在所有工況條件以及催化劑幾何性質均相同的條件下，三葉草型催化劑內部的最高溫度及平均溫度比相對應的圓柱型催化劑一般低0.1～0.9 ℃，而且其催化劑內部溫度波動范圍也要小于常規(guī)圓柱型催化劑。

表2所對應的相關不同工藝條件其催化劑內部的溫度分布比較相似，因此將15號計算對應的三維溫度場分布做二維軸向及徑向切面圖，結果如圖6所示。

表2 熱點分布均勻的三葉草型催化劑表面溫度370 ℃及不同工況下的內部溫度

(a) 整體

(b) 軸向

(c) 徑向圖6 序號15條件三葉草型催化劑內部溫度分布切片圖

對于熱點分布均勻的等溫反應條件，由于三葉草型催化劑的截面較常規(guī)圓柱型催化劑復雜，因此其內部的溫度分布形式也與常規(guī)的圓柱型催化劑相差較大，尤其是體現(xiàn)出了較為復雜的赤道角向信息[13]。分析圖6可知，溫度分布的等勢面在三葉草型催化劑截面每個葉上的徑向方向呈近似貝殼型的分布，沿著大圓圓心與小圓圓心連線方向其溫度梯度變化最小，而沿著大圓圓心與小圓共切點方向其溫度梯度變化最大，由于催化劑體系具有軸對稱的特征，其最高溫度與圓柱型催化劑相同，仍在催化劑的中心軸線的中點位置上出現(xiàn)。而在軸向剖面上其溫度等勢面呈現(xiàn)近似圓角矩形分布，在外表面溫度為370 ℃的催化劑內部溫度逐漸增加，其中徑向方向溫度梯度大于軸向方向。也可以認為雖然在理論上催化劑內部并非等溫環(huán)境，但其內部的反應環(huán)境與等溫環(huán)境比較類似。但當催化劑半徑及長度較大時，其三葉草型催化劑中心最高溫度與表面溫度相差較大，在該催化劑內部的反應相對于等溫反應偏離較大。因此，僅就熱效應而言，制備三葉草型催化劑時其半徑及長度一般也不宜過大。

3 結 論

通過本文的綜合分析可以發(fā)現(xiàn)，即使在理想的宏觀等溫反應條件下，圓柱型及三葉草型催化劑內部仍然是具有非恒定溫度場分布的非等溫反應區(qū)域。在常規(guī)加氫裂化反應中，較高的反應熱、較高的處理空速、較大的原料油密度以及較大的催化劑半徑及催化劑軸向長度均導致較高的催化劑內部溫度。其中反應熱和催化劑半徑對于催化劑內部溫度的影響更顯著。同時，在其他條件相同時，三葉草型催化劑內部的最高溫度、平均溫度均低于圓柱型催化劑，同時其擁有比圓柱型催化劑相對更均勻的溫度分布。一般加氫催化反應過程中，催化劑內部的熱點分布對于催化劑內部的溫度分布影響比較顯著。因此，確立比較接近工業(yè)催化劑實際環(huán)境的熱源函數(shù)對于催化劑內部溫度的進一步分析十分必要。

[1] 尹兆林.煉油企業(yè)全廠用能分析及節(jié)能優(yōu)化[J].石油煉制與化工，2012，43(10):86-91.

[2] 董兆海, 袁永新,王明傳.加氫裂化裝置能耗及節(jié)能分析[J].齊魯石油化工, 2011,39(2):87-91.

[3] 劉小波，毛羽，王娟，等.基于多孔介質加氫裂化反應器多相流數(shù)值模擬[J].石油學報(石油加工)，2012,28(2):260-267.

[4] 鄭忠，何臘梅.紅外測溫技術及在鋼鐵生產(chǎn)中的應用[J].工業(yè)加熱，2005, 34(3):25-29.

[5] 廖玉麟.數(shù)學物理方程[M].上海:華東理工大學出版社，1995：12-13.

[6] 龍軍, 邵潛, 賀振富, 等.規(guī)整結構催化劑及反應器研究進展[J].化工進展, 2004,23(9):925-931.

[7] 梁文杰，闕國和.石油化學[M].東營:中國石油大學出版社，2011:356.

[8] 吳建民， 張海濤， 應衛(wèi)勇， 等.鈷基催化劑固定床有效導熱系數(shù)[J].過程工程學報，2010,10(1):29-34.

[9] 劉欣.無網(wǎng)格方法[M].北京：科學出版社，2011：25-29.

[10]LiuGR，GuYT.無網(wǎng)格法理論及程序設計[M].王建明，周學軍，譯.濟南:山東大學出版社，2007:24.

[11] 劉維.實戰(zhàn)Matlab之并行程序設計[M].北京：北京航空航天大學出版社，2012:135-137.

[12] 王曉丹，吳崇明.基于MATLAB的系統(tǒng)分析與設計：圖像處理[M].西安:西安電子科技大學出版社，2000:7.

[13] 莫爾比代利M，加夫里迪斯A，瓦爾馬A.催化劑設計[M].北京:化學工業(yè)出版社，2004:128.

(編輯 宋官龍)

Numerical Simulation about Bulk Temperature Distribution of Clover-Type Catalyst

Zhao Bo1, Wang Kuo2， Sui Chengyu1

(1.CollegeofScience,LiaoningShihuaUniversity,FushunLiaoning113001,China;2.FushunResearchInstituteofPetroleumandPetrochemical,FushunLiaoning113001,China)

Based on three-dimensional somatic environment of the real clover type hydrocracking catalyst as the computational entity and the simulated industrial operating temperature as the boundary condition. The Fourier heat transfer equation is solved by the meshless numerical method. Then the influence of external temperature fluctuation on the internal temperature distribution of the catalyst was analyzed by using the calculated results. The analysis results show that the actual reactions in the catalyst is not isothermal reaction. At the same time, the fluctuation of the reaction temperature and the hot spot inside the catalyst have a certain influence on the temperature distribution for the catalyst cluster in the hydrocracking process. Catalyst bulk phase inner average temperature is influenced by heat of reaction, the catalyst particle size, material density, reaction space velocity and catalyst inner hotspot distribution.

Hydrocracking; Heat transfer; Meshless; Numerical simulation; Distributed computing

1672-6952(2017)02-0010-07

2016-09-09

2016-09-29

中國石油化工股份有限公司項目(011308)。

趙波(1980-)，女，碩士，講師，從事物理實驗和數(shù)值模擬研究；E-mail：381203658@qq.com。

TQ012

A

10.3969/j.issn.1672-6952.2017.02.003

投稿網(wǎng)址：http://journal.lnpu.edu.cn